Predictions of Tropical Forest Biomass and Biomass Growth Based on Stand Height or Canopy Area Are Improved by Landsat-scale Phenology across Puerto Rico and the U.s. Virgin Islands

Remotely-sensed estimates of forest biomass are usually based on various measurements of canopy height, area, volume or texture, as derived from LiDAR, radar or fine spatial resolution imagery. These measurements are then calibrated to estimates of stand biomass that are primarily based on tree stem diameters. Although humid tropical forest seasonality can have low amplitudes compared with temperate regions, seasonal variations in growth-related factors like temperature, humidity, rainfall, wind speed and day length affect both tropical forest deciduousness and tree height-diameter relationships. Consequently, seasonal patterns in spectral measures of canopy greenness derived from satellite imagery should be related to tree height-diameter relationships and hence to estimates of forest biomass or biomass growth that are based on stand height or canopy area. In this study, we tested whether satellite image-based measures of tropical forest phenology, as estimated by indices of seasonal patterns in canopy greenness constructed from Landsat satellite images, can explain the variability in forest deciduousness, forest biomass and net biomass growth after already accounting for stand height or canopy area. Models of forest biomass that added phenology variables to structural variables similar to those that can be estimated by LiDAR or very high-spatial resolution imagery, like canopy height and crown area, explained up to 12% more variation in biomass. Adding phenology to structural variables explained up to 25% more variation in Net Biomass Growth (NBG). In all of the models, phenology contributed more as interaction terms than as single-effect terms. In addition, models run on only fully-forested plots performed better than models that included partially-forested plots. For forest NBG, the models produced better results when only those plots with a positive growth, from Inventory Cycle 1 to Inventory Cycle 2, were analyzed, as compared to models that included all plots.


Introduction
Forest carbon offset programs, like the United Nations program for Reducing Emissions from Deforestation and forest Degradation (REDD), which is aimed at tropical forests, have the potential to mitigate human-caused increases in atmospheric greenhouse gases while helping to conserve the high biodiversity of tropical forests by providing financial incentives to reduce atmospheric greenhouse gas emissions from forests.The REDD program requires inventory and monitoring of tropical forest carbon stocks and their dynamics [1], the aboveground live components of which are generally estimated from tree biomass and changes in tree biomass.In the context of REDD, our research targets the aboveground live biomass carbon component of Equation 2.3 of the IPCC's 2006 Guidelines for National Greenhouse Gas Inventories [2]: where ∆C LU is the carbon stock changes for a land-use stratum, and the subscripts denote the following carbon pools over the same time period: AB: aboveground biomass BB: belowground biomass DW: deadwood LI: litter SO: soils HWP: harvested wood products.
Forest aboveground live biomass is most commonly estimated by summing tree biomass, estimated from tree stem diameter, as measured in forest inventory plots, with allometric equations [3], though tree height or wood density may also be included in calculations.Remote sensing estimates of forest biomass are usually derived from LiDAR or very fine spatial resolution imagery (e.g., with pixel sizes of 1 m or less) and are based on calibrating measures of remotely-sensed canopy height, area, volume or texture, to field-based estimates of forest biomass (e.g., Couteron et [4][5][6][7]).Rainfall affects not only tropical forest phenology [8,9], but also tree height-diameter relationships [3,10].Consequently, seasonal patterns in spectral measures of canopy greenness derived from satellite imagery, which we refer to here as phenology, should be related to tree height-diameter relationships and, hence, to estimates of forest biomass or biomass growth.As a result, phenology data have the potential to improve forest biomass estimates based on canopy characteristics, such as those that might be quantified with LiDAR or very fine spatial resolution imagery.
Detailed phenology products can be obtained from near-surface remote sensing [11][12][13], but such intensive monitoring of vegetation phenology is probably not practical for REDD.Because of persistent cloud cover in tropical areas, satellites with high temporal resolution, having revisit times of one to two days, such as the Moderate Resolution Imaging Spectroradiometer (MODIS), Advanced Very High Resolution Radiometer (AVHRR) or the Greenhouse Gases Observing Satellite (GOSAT), are more likely than satellites with less frequent coverage to collect enough cloud-free images to track annual changes in forest greenness.Vegetation phenology studies based on such short revisit time imagery use indices like the Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Chlorophyll/Carotenoid Index (CCI) or Sun-Induced Chlorophyll Fluorescence (SIF) [14][15][16][17][18][19].Although these products have a high temporal frequency, they have a coarse spatial resolution ranging from 250 m to several kilometers.This range of spatial scales is too coarse to characterize the spatial variability in forest attributes across complex forest landscapes, where biophysical attributes and forest types can change over much shorter distances.Moreover, such phenology products do not precisely reflect the phenology of most field plots, because their pixel size is larger than most field plots, which are as small as 0.067 ha for the U.S. Forestry Inventory and Analysis Program.Furthermore, large pixels are much more likely to reflect the phenology of both forest and non-forest land cover.
Phenology indices from satellite imagery with finer spatial resolution, such as Landsat imagery, are a potential alternative to frequent, but coarsely-scaled imagery or near-surface methods for gauging forest phenology.We already know that that multi-season Landsat satellite imagery is useful for distinguishing among tropical forest types with different degrees of deciduousness, both visually and in automated classifications [20,21].Phenology data at Landsat resolution should, therefore, provide a product that minimizes the uncertainty related to mixtures in land cover and forest types and is closer to the size of forest field plots.The main challenge with Landsat is that its 16-day revisit time reduces the availability of cloud-free images, calling into question whether phenology data derived from Landsat would have too much cloud contamination to be useful.Furthermore, investigating this question is difficult as it requires, at a minimum, a large investment in cloud and shadow masking.Probably, for this reason, no previous studies have investigated whether Landsat-derived phenology in a persistently cloudy tropical landscape can improve remote sensing estimates of biomass or biomass growth.
In this study, we test whether indices related to greenness seasonality, as estimated from repaired Landsat imagery (i.e., missing data due to clouds, shadow and Scan Line Corrector (SCL) off stripping, have been repaired) will help explain variability not only in the relative importance of deciduous species, but also in aboveground live biomass and biomass growth, even after accounting for stand height and canopy area.If phenology can explain additional variability in forest biomass or growth, then the potential exists for forest phenology measures to reduce uncertainty in remotely-sensed estimates of biomass or growth from LiDAR or fine spatial resolution imagery.To test this hypothesis, we used a set of new algorithms [22,23] that detect and mask clouds and shadow, repair images and reconstruct seasonal Landsat time series data and detect tropical forest phenology metrics at the spatial resolution of Landsat data, which is 30 m.We then evaluated the extent to which phenology information contributed to explaining the variation in forest biomass and net biomass growth after accounting for forest stand height or canopy area along a gradient from tropical dry to moist and wet forests in Puerto Rico and the U.S. Virgin Islands (PRVI).We also evaluated whether the phenology indices are related to the relative importance of deciduous species.PRVI is an ideal setting to test this hypothesis since its forests range from dry to wet, and rainfall and temperature seasonality affect forest species composition, including the relative importance of deciduous species [24].

Study Area
Puerto Rico and the United States Virgin Islands are biologically-diverse tropical islands in the Caribbean Sea.We found this to be a good site to test our hypotheses because of the vegetation biophysical diversity that results from substantial gradients in soil types, temperature and rainfall, at an observational scale of the Forest Inventory and Analysis (FIA) plots that corresponds well with Landsat pixels.The mean annual precipitation ranges from <700 mm in the southwest to >4500 mm in the northeast [25].The sites are rich in forest type and tree species diversity due to the diverse climatic zones, previous land use types and soils that have formed on alluvial, sedimentary, volcanic, limestone or serpentine substrates [26].The most extensive climatic zone has broadleaf evergreen forest that on karst substrate forms a mosaic of evergreen forest and forest that is semi-deciduous during the dry season [26].Other forest types include dry semi-deciduous forest and cloud forest (at the highest elevations above the cloud condensation level).

Methods
The methodological framework for this work is shown in Figure 1.We used field data collected by the Forest Inventory and Analysis (FIA) Program (Section 2.2.1) and remote sensing data from the Landsat ETM+ sensor.The FIA Program is conducted by the U.S. Department of Agriculture Forest Service and utilizes a nationally-standardized sampling design for monitoring the forest of the U.S., including PRVI.Systematically installed field plots across the lands of PRVI are revisited every five years in what is referred to as one cycle of Caribbean FIA data.Each permanent field plot represents a cluster of four 7.3-m radius subplots (total sampled area = 0.067 ha), where a standard set of variables are collected to describe land use, stand structure, tree species composition and regeneration on all forested sites during each cycle.Trees with a minimum diameter at breast height (dbh) of 12.7 cm are surveyed for each subplot.Each subplot includes a micro-plot with a 2.1-m radius where trees with a dbh from 2.54-12.7 cm are surveyed.For more information regarding the FIA Program and plot design, see Bechtold and Patterson [27].Landsat ETM+ images were used as they are the only source of Landsat data for this area over the 2006-2008 study period.PRVI lies in one of the tropical areas that have a mean cloud fraction in Landsat ETM+ acquisition exceeding 60% [28].Such cloud contamination reduces image usability for land surface studies; thus, we needed to reconstruct cloud-free Landsat images.The reconstruction involves two steps: automatic detection of clouds and shadows, and interpolation of missing data caused by clouds, shadows and the lines or missing data caused by the failure of the Landsat 7 ETM+ Scan Line Corrector.In both steps, our study applied cutting-edge algorithms to achieve high accuracy.The reconstructed images were then used for detecting phenology (Section 2.2.2).

Field Data Collection and Processing
We used field data collected during two FIA cycles: inventory years 2001-2004 (t1) and inventory years 2006-2010 (t2).Height, deciduousness proportions and aboveground biomass data were available for both cycles, but crown area was available for t1 data only.In this paper, "all plots" refers to those FIA plots with one to four forested subplots, while "100% forest plots" means those plots in which all four subplots were forested.Aboveground biomass is estimated by the FIA Program in PRVI by using locally-developed allometric equations for forest life zones sensu [3,[29][30][31][32][33]. Landsat ETM+ images were used as they are the only source of Landsat data for this area over the 2006-2008 study period.PRVI lies in one of the tropical areas that have a mean cloud fraction in Landsat ETM+ acquisition exceeding 60% [28].Such cloud contamination reduces image usability for land surface studies; thus, we needed to reconstruct cloud-free Landsat images.The reconstruction involves two steps: automatic detection of clouds and shadows, and interpolation of missing data caused by clouds, shadows and the lines or missing data caused by the failure of the Landsat 7 ETM+ Scan Line Corrector.In both steps, our study applied cutting-edge algorithms to achieve high accuracy.The reconstructed images were then used for detecting phenology (Section 2.2.2).

Field Data Collection and Processing
We used field data collected during two FIA cycles: inventory years 2001-2004 (t 1 ) and inventory years 2006-2010 (t 2 ).Height, deciduousness proportions and aboveground biomass data were available for both cycles, but crown area was available for t 1 data only.In this paper, "all plots" refers to those FIA plots with one to four forested subplots, while "100% forest plots" means those plots in which all four subplots were forested.Aboveground biomass is estimated by the FIA Program in PRVI by using locally-developed allometric equations for forest life zones sensu [3,[29][30][31][32][33]. Species-specific aboveground biomass equations are exceptionally available and used for species such as Prestoea montana and Bucida buceras [3,30,33].Tables showing the equations used by FIA to predict aboveground biomass in PRVI are available in Brandeis et al. [34] and Brandeis and Oswalt [35].Net biomass growth was calculated as follows: where all biomass values are given in Mg•ha −1 of live biomass, dry weight and: t 1 = Julian day of the first survey in fractional y t 2 = Julian day of the second survey in fractional y t m = Julian day of the midpoint between the first and second surveys in fractional y t 2 − t 1 = Interval between the first and second survey dates in fractional y NBG = Net Biomass Growth at the plot level, in Mg•ha −1 •y −1 , from t 1 to t 2 SAB tx = Biomass of survivor tree i estimated for t x , as indicated IAB tx = Biomass of ingrowth tree i estimated for t x , as indicated MAB tx = Biomass of mortality tree i estimated for t x , as indicated S = The number of trees that survived from t 1 to t 2 I = The number of trees in micro-plots that reached a dbh of 12.7 cm in t 2 M = The number of trees that died in the interval from t 1 to t 2 Total aboveground tree biomass and net biomass growth for each plot were estimated using the total plot area, including both forest and non-forest plot area.This approach relies on the idea that in the 3 × 3 window from which we extracted Landsat-derived variables, non-forest areas would reduce both the image variables and the forest variables similarly and in proportion with non-forest area.For field measured independent variables, we chose to use only those that have been proven to be readily derivable from remote sensing data, such as LiDAR, radar or high-spatial resolution multispectral images.We estimated crown areas from crown diameters (which are measured in two perpendicular directions), assuming that crowns had an elliptical shape.Three main plot height variables were calculated: the maximum tree height (H max ), the mean height of all live dominant and co-dominant trees (H domcod ) and crown area weighted height (H cawgt ).We weighted height by crown area rather than by stem basal area (as is done when estimating Lorey's height) assuming that measures of crown area or volume might be derived from LiDAR or high-spatial resolution remote sensing, but not stem diameters, which are typically used to calibrate the remote sensing data to biomass as estimated with field-measured stem diameter.One study that estimated stand-level stem diameter at breast height with LiDAR-derived variables [36] found that they only explained about half of its variation.We dropped three forest plots with strong negative net biomass growth due to very high mortality, or because for one large tree, height at t 2 was much shorter than height at t 1 , even though no damage was recorded.

Phenology Detection
We reconstructed a time series dataset with Landsat 7 ETM+ images covering PRVI.We only selected images with cloud coverage lower than 50% because the useful information in images with higher cloud coverage is very limited.An annual phenology product ideally requires one year's worth of time series images, which is impossible to get in a cloudy tropical area.To address this, our approach uses images from as few consecutive years as possible, ordered by Julian day of the year (DOY) and treats them as a one-year time series [23].Our approach tests the idea that, even though forest phenology changes from year to year, there are enough similarities among years to form a pattern that is indicative of forest characteristics related to biomass and biomass growth.We used images from 2006-2008, corresponding to the time period when most of the latest field data were collected.Four Landsat scenes (Path 4, Rows 47-48 and Path 5, Rows 47-48) are required to cover PRVI (excluding Mona Island).After filtering for cloud cover, a minimum of 17 and maximum of 28 Landsat images were available from each scene to reconstruct the time series (Table 1).In this study, we applied a newly-developed algorithm to screen clouds and cloud shadows in all Landsat images automatically [22].This new algorithm makes better use of the temporal information in the time series and can obtain higher accuracy than other algorithms [37][38][39][40][41].After the detection of clouds and shadows, we used an advanced spatiotemporal interpolation method, named the nearest similar pixel interpolator (NSPI), to estimate the values of cloud and shadow pixels [42] and un-scanned pixels within SLC-off gaps [43].Chen et al. [44] developed the weighted cross-correlogram spectral matching-phenology (CCSM-P) method that detects phenology from vegetation index time series of multiple years.Zhu et al. [23] modified CCSM-P so that it can also detect phenology for a single year time series.From the reconstructed images, we computed a time series of the Normalized Difference Vegetation Index (NDVI).We believe NDVI is a simple index suitable for a first product like ours; thus, the necessity to use alternative and even more complicated indices depends on the output of the work presented here.A time series smoothing algorithm [45] based on the Savitzky-Golay filter [46] was used to further remove noise in the NDVI time series from which we identified 5 main variables: minimum NDVI (ndvi_min), maximum NDVI (ndvi_max), NDVI range (ndvi_range), mean NDVI (ndvi_mean) and NDVI standard deviation (ndvi_std).The modified CCSM-P algorithm [23] was adopted for the smoothed NDVI time series data to detect 5 important phenology variables: the DOY with the highest NDVI value (peak), DOY with the lowest NDVI (lowest), the DOY when growth starts (greenup), the DOY when growth starts to slow down (browndown) and the length of time between browndown and greenup time, dry season length (dsl).Using the FIA plot center coordinates, each plot was assigned phenology variables by averaging the values of its surrounding 3 × 3 pixel window, thus approximating the main plot's dimensions.

Statistical Analysis
To test the potential contribution of phenology data in explaining forest biomass, we used data from t 2 to develop a simple linear model (Equation ( 2)) relating biomass to canopy height and used the model's residuals to determine the annual phenology pattern of high biomass (upper quartile) versus low biomass (lower quartile) plots.We then developed another linear model (Equation ( 3)) relating biomass to canopy height and forest deciduousness (proportional area covered by deciduous tree species, including deciduous, semi-deciduous or facultatively deciduous trees).
where AGB is aboveground biomass in mega-grams per hectare, Ht is crown area weighted height in meters, Dec is the proportion of plot basal area covered by deciduous trees and a, b, c, d and e are constants.We used crown area weighted height because it produced the best model (highest variation in biomass explained) among the three height variables mentioned in Section 2.2.1.
From the two models, we compared the annual pattern of phenology from the residuals of the second model against the pattern from the first model (Figure 2).This figure shows that lower biomass (lower residuals) plots have a larger annual range/variability in NDVI, but the trends become very similar after we use the deciduousness proportions, thereby demonstrating that phenology data might explain some of the variation in forest physiognomy.The statistical contribution of phenology data in modeling forest deciduousness was then tested by developing a multiple linear regression model (Equation ( 4)) that estimated forest deciduousness proportion from both DOY and NDVI phenology variables.For forest biomass and net biomass growth, we initially tested the relationship between aboveground biomass and canopy height or crown area for a single inventory period (tn) using simple ordinary least squares regression.We then implemented forward stepwise multiple linear regression analysis to determine how much additional variation is explained by adding phenology variables.To also determine the effect of forest type proportions, we first did the analysis on all FIA plots and then only on the 100% forest plots and compared the two results.Since we had biomass data for the t1 to t2 (Δt) growth period and with the argument that we can get specialized remote sensing data to estimate canopy height or crown area at tn and phenology data to estimate growth patterns, we tested the ability of height at tn and delta height (ΔHt) plus phenology variables to estimate NBG.We could have also done the same for crown area, but we did not have this variable for t2 data.For the multiple linear regression processes, each final model was run using only those variables that were significant at the α = 0.05 level, and we recorded the amount of variation explained as denoted by the adjusted R 2 value.
We also performed a quality check on the gap/cloud/shadow repairing process.For each plot, we checked how many of the 9 pixels in the 3 × 3 window had original image values and how many had filled (repaired) values for each of the dates.We then calculated each date's "quality" as the number of pixels with original values divided by 9, thus to get the overall quality for each plot, we averaged all of the dates' quality values for that respective plot.For any model, the residuals did not show any trend with the overall quality, suggesting that our results are not heavily influenced by the final image quality.Thus, we concluded that we repaired the SLC-off and cloud/shadow gaps well and, hence, addressed the problem of cloud cover that often limits the use of medium temporal resolution imagery in studying the phenology of tropical forests [1,47].

Phenology
As Figure 3 and Table 2 show, the annual range of NDVI is very low (mean of 0.15).On average, the forests have the lowest greenness (mean NDVI of 0.64) in April-May and peak greenness (mean For forest biomass and net biomass growth, we initially tested the relationship between aboveground biomass and canopy height or crown area for a single inventory period (t n ) using simple ordinary least squares regression.We then implemented forward stepwise multiple linear regression analysis to determine how much additional variation is explained by adding phenology variables.To also determine the effect of forest type proportions, we first did the analysis on all FIA plots and then only on the 100% forest plots and compared the two results.Since we had biomass data for the t 1 to t 2 (∆t) growth period and with the argument that we can get specialized remote sensing data to estimate canopy height or crown area at t n and phenology data to estimate growth patterns, we tested the ability of height at t n and delta height (∆Ht) plus phenology variables to estimate NBG.We could have also done the same for crown area, but we did not have this variable for t 2 data.For the multiple linear regression processes, each final model was run using only those variables that were significant at the α = 0.05 level, and we recorded the amount of variation explained as denoted by the adjusted R 2 value.
We also performed a quality check on the gap/cloud/shadow repairing process.For each plot, we checked how many of the 9 pixels in the 3 × 3 window had original image values and how many had filled (repaired) values for each of the dates.We then calculated each date's "quality" as the number of pixels with original values divided by 9, thus to get the overall quality for each plot, we averaged all of the dates' quality values for that respective plot.For any model, the residuals did not show any trend with the overall quality, suggesting that our results are not heavily influenced by the final image quality.Thus, we concluded that we repaired the SLC-off and cloud/shadow gaps well and, hence, addressed the problem of cloud cover that often limits the use of medium temporal resolution imagery in studying the phenology of tropical forests [1,47].

Phenology
As Figure 3 and Table 2 show, the annual range of NDVI is very low (mean of 0.15).On average, the forests have the lowest greenness (mean NDVI of 0.64) in April-May and peak greenness (mean NDVI of 0.80) in November-December.Browning starts around February-March, while greening up starts around June-July, giving a dry season length of 4-5 months.Note that the dry season length is relative to the PRVI conditions, thus not necessarily dry, but rather, a period when greenness is below the average.The spatial pattern of phenology is shown in Figures 4 and 5.
Remote Sens. 2017, 9, 123 8 of 18 NDVI of 0.80) in November-December.Browning starts around February-March, while greening up starts around June-July, giving a dry season length of 4-5 months.Note that the dry season length is relative to the PRVI conditions, thus not necessarily dry, but rather, a period when greenness is below the average.The spatial pattern of phenology is shown in Figures 4 and 5.   NDVI of 0.80) in November-December.Browning starts around February-March, while greening up starts around June-July, giving a dry season length of 4-5 months.Note that the dry season length is relative to the PRVI conditions, thus not necessarily dry, but rather, a period when greenness is below the average.The spatial pattern of phenology is shown in Figures 4 and 5.

Forest Deciduousness
Directly relating forest deciduousness to phenology variables showed weak, but significant (p < 0.001) correlations (Figure 6).When all plots were considered, greenup date, minimum NDVI and the interaction of greenup date and minimum NDVI explained 26% of the variation in forest deciduousness.When only 100% of the forest plots were used in the model, the amount of variation explained was 36%, and the significant variables were greenup date, lowest date, standard deviation of NDVI and the interaction of browndown date and standard deviation of NDVI.

Forest Deciduousness
Directly relating forest deciduousness to phenology variables showed weak, but significant (p < 0.001) correlations (Figure 6).When all plots were considered, greenup date, minimum NDVI and the interaction of greenup date and minimum NDVI explained 26% of the variation in forest deciduousness.When only 100% of the forest plots were used in the model, the amount of variation explained was 36%, and the significant variables were greenup date, lowest date, standard deviation of NDVI and the interaction of browndown date and standard deviation of NDVI.

Forest Biomass
Table 3 summarizes the contribution of phenology in explaining the variation in biomass.Phenology added a reasonable amount of variation explained by the models, more as interaction terms than as single effect terms.The height variables tested performed similarly in most situations.Crown area-based models performed slightly better than height-based models, but models with both crown area and height performed much better than models that used one of the variables alone.Both NDVI and DOY variables were useful, and the combination depended on the type of model and the group of forest plots.Overall, models for 100% forest plots had better results than models for all plots combined, probably due to the lower variability in the former.Table 3 shows the details of the models for each group of plots, and Figure 7 shows the best models for each of the groups.

Forest Biomass
Table 3 summarizes the contribution of phenology in explaining the variation in biomass.Phenology added a reasonable amount of variation explained by the models, more as interaction terms than as single effect terms.The height variables tested performed similarly in most situations.Crown area-based models performed slightly better than height-based models, but models with both crown area and height performed much better than models that used one of the variables alone.Both NDVI and DOY variables were useful, and the combination depended on the type of model and the group of forest plots.Overall, models for 100% forest plots had better results than models for all plots combined, probably due to the lower variability in the former.Table 3 shows the details of the models for each group of plots, and Figure 7 shows the best models for each of the groups.

Forest Biomass
Table 3 summarizes the contribution of phenology in explaining the variation in biomass.Phenology added a reasonable amount of variation explained by the models, more as interaction terms than as single effect terms.The height variables tested performed similarly in most situations.Crown area-based models performed slightly better than height-based models, but models with both crown area and height performed much better than models that used one of the variables alone.Both NDVI and DOY variables were useful, and the combination depended on the type of model and the group of forest plots.Overall, models for 100% forest plots had better results than models for all plots combined, probably due to the lower variability in the former.Table 3 shows the details of the models for each group of plots, and Figure 7 shows the best models for each of the groups.

Forest Net Biomass Growth
As Table 4 and Figure 8 show, phenology variables were important in the NBG models mostly in the same manner as they were for forest biomass models.Furthermore, the models for 100% forest plots were generally better than the models for all plots combined.The models produced better results when only forests with positive growth were considered, thus suggesting a higher complexity in mortality drivers in the negative growth plots compared to growth drivers in the positive growth plots.

Forest Net Biomass Growth
As Table 4 and Figure 8 show, phenology variables were important in the NBG models mostly in the same manner as they were for forest biomass models.Furthermore, the models for 100% forest plots were generally better than the models for all plots combined.The models produced better results when only forests with positive growth were considered, thus suggesting a higher complexity in mortality drivers in the negative growth plots compared to growth drivers in the positive growth plots.

Discussion
This study has demonstrated how phenology from moderate to fine spatial resolution imagery helps quantify plot-level forest biomass, growth and deciduousness in highly cloudy forest environments.We have revealed that despite the heterogeneity and complexity of tropical forest ecosystems, phenology data significantly improve the estimation of both t n and ∆t biophysical attributes.The ability of phenology data to quantify the growth pattern of vegetation complements the ability of structural variables to estimate the change in the structure of vegetation over time and, thus, is a good indicator of its function, such as carbon sequestration or a reflection of its response to growth-controlling factors, such as climate change.The complementary contribution of phenology is reflected by the significance of mostly interactive variables as compared to single effect variables in the more complex models.The better performance of models when only 100% forested plots were considered strongly suggests that pixels with mixed land cover, which would be more common in coarse spatial resolution imagery like MODIS, adding noise to the relationships between phenology and the forest attributes that we tested.
Phenology variables like NDVI range or growing season length are expected to clearly distinguish the deciduous proportions of forest plots; thus, we suspect that the mixture of old dominant and new sub-canopy vegetation in PRVI complicated phenology patterns, hence the modest results in our forest deciduousness models.The two most common and important species, like Guarea guidonia and Spathodea campanulata, even have a bimodal flowering/leaf fall pattern [34,[48][49][50].Other species, like Tabebuia heterophylla, vary in their deciduousness, depending on climate.These complexities probably contributed to the modest relationship between the phenology indices and proportion of basal area composed of deciduous species.This complexity could also be the reason why extreme-end phenology variables like ndvi_min, ndvi_max and peak DOY were mostly significant in the biomass and NBG models, compared to range-based variables.Depending on the amount of photosynthetically-active radiation available in the form of diffuse light, the understory may experience different seasonality in irradiation and, hence, phenology compared to the canopy [51].Furthermore, the high diversity of plant species in PRVI, combined with differences in their phenological patterns, may lead to some unexpected patterns in seasonality at the community and landscape level [52], posing a challenge in detecting phenology from moderate to coarse spatial resolution remote sensing data.For future studies, improvements may come from incorporating stand age as an explanatory variable or using alternative or more vegetation indices in the models.Fan et al. [53], for example, used Landsat 8 time series of NDVI, EVI, the Atmospherically-Resistant Vegetation Index (ARVI), Normalized Difference Moisture Index (NDMI) and Tasseled Cap Greenness (TCG) to identify foliation (leaf flushing) and defoliation (leaf-off) days to help distinguish rubber trees from natural forests and croplands with up to 96% (kappa = 0.92) accuracy in Tropical South West China.For a first cut analysis, we chose NDVI as a commonly-used measure of greenness.Now that we obtained promising results, it will be interesting to test more robust and structurally sensitive indices in future work.A longer time series dataset may also improve the phenology variables if the associated land cover dynamics are also carefully controlled for.
It is important to note that net biomass growth models can be much better when more variables that account for climate variability, nutrients or disturbance/mortality are used.Reich [54] modeled annual and daily net primary productivity (NPP) from canopy traits (Leaf Area Index (LAI) and percentage nitrogen (N)) and climate variables and also reported the importance of interaction terms in improving the amount of variation explained.A model that included LAI, N, summer rainfall, mean maximum summer temperature and growing season length explained 92% of the variation in annual NPP, an additional 16% to a model that used LAI and N alone.However, the main challenge of including these variables, such as percentage nitrogen, is that they are not generally derivable from multispectral remote sensing data, which is what is normally available for a region.Moreover, our main objective was to test the potential contributions of Landsat level phenology, given the wide availability of Landsat.
Considering the accuracy of the phenology characterization, in terms of the resulting times, for MODIS-based EVI time series data, Zhang et al. [55] in North American temperate forests concluded that phenology can be detected with an accuracy of less than three days for a temporal resolution of less than 16 days and that the error increases with the number of missing dates, especially if these missing days are around the onset of phenological transition dates.Compared with phenology characterizations where a denser time series of images is available, this study necessarily uses fewer images, as few as 17 images for the Path 5, Row 47 scene.Although using fewer images may have reduced accuracy, the time series reconstruction algorithm that we used is designed to correct this problem by using a time-based function rather than an observation points-based interpolator [23,45].However, our study focused more on how much variability in forest biomass, net biomass growth and deciduousness was explained by phenology rather than the absolute accuracy of the phenology indices.

Conclusions
This study demonstrated the possibility of using moderate temporal resolution, but high spatial resolution time series images to generate phenology variables that are important for estimating forest aboveground live biomass and net biomass growth across a tropical landscape of high diversity and complexity.Monitoring phenology in tropical forests has been a great challenge using the traditional medium to coarse temporal resolution multispectral imagery due to residual cloud and aerosol contamination.However, as we have demonstrated in this paper, the experience gained in using the traditional high-frequency, coarse spatial resolution imagery is vital for creating cloud and gap-filled medium-to high-spatial resolution imagery for monitoring phenology at local scales.
Remote Sens. 2017, 9, 123 4 of 18 variables are collected to describe land use, stand structure, tree species composition and regeneration on all forested sites during each cycle.Trees with a minimum diameter at breast height (dbh) of 12.7 cm are surveyed for each subplot.Each subplot includes a micro-plot with a 2.1-m radius where trees with a dbh from 2.54-12.7 cm are surveyed.For more information regarding the FIA Program and plot design, see Bechtold and Patterson [27].

Figure 1 .
Figure 1.The data and methodological steps followed in phenology detection and final forest deciduousness, forest biomass and forest net biomass growth modeling.

Figure 1 .
Figure 1.The data and methodological steps followed in phenology detection and final forest deciduousness, forest biomass and forest net biomass growth modeling.

Forest Deciduousness = 4 . 4 ) 4 )Figure 2 .
Figure 2. Variability in phenology trends according to biomass levels and model residuals.Left: How mean monthly NDVI variability is related to the residuals of a model that uses height only.Right: How mean monthly NDVI variability is related to the residuals of a model that used height and forest deciduousness variables.

Figure 2 .
Figure 2. Variability in phenology trends according to biomass levels and model residuals.Left: How mean monthly NDVI variability is related to the residuals of a model that uses height only.Right: How mean monthly NDVI variability is related to the residuals of a model that used height and forest deciduousness variables.

Figure 3 .
Figure 3. Annual temporal variability in mean plot-level NDVI for in Puerto Rico and the U.S. Virgin Islands (PRVI) Forest Inventory and Analysis (FIA) field plots.

Figure 4 .
Figure 4. Spatial variability of phenology expressed in NDVI variables.

Figure 3 .
Figure 3. Annual temporal variability in mean plot-level NDVI for in Puerto Rico and the U.S. Virgin Islands (PRVI) Forest Inventory and Analysis (FIA) field plots.

Figure 3 .
Figure 3. Annual temporal variability in mean plot-level NDVI for in Puerto Rico and the U.S. Virgin Islands (PRVI) Forest Inventory and Analysis (FIA) field plots.

Figure 4 .
Figure 4. Spatial variability of phenology expressed in NDVI variables.

Figure 4 .
Figure 4. Spatial variability of phenology expressed in NDVI variables.

Figure 5 .
Figure 5. Spatial variability of phenology expressed in Julian day of the year variables.

Figure 5 .
Figure 5. Spatial variability of phenology expressed in Julian day of the year variables.

Figure 6 .
Figure 6.Modeling forest deciduousness at t2: proportion of deciduous cover as a function of phenology.

Figure 7 .
Figure 7. Results of the best models for forest biomass.For each row, the left two columns are for a model run on all plots, and the right two columns are for a model run on 100% forest plots only.

Figure 6 .
Figure 6.Modeling forest deciduousness at t 2 : proportion of deciduous cover as a function of phenology.

Figure 6 .
Figure 6.Modeling forest deciduousness at t2: proportion of deciduous cover as a function of phenology.

Figure 7 .
Figure 7. Results of the best models for forest biomass.For each row, the left two columns are for a model run on all plots, and the right two columns are for a model run on 100% forest plots only.

Figure 7 .
Figure 7. Results of the best models for forest biomass.For each row, the left two columns are for a model run on all plots, and the right two columns are for a model run on 100% forest plots only.

Figure 8 .
Figure 8. Results of the best models for forest net biomass growth (NBG).Top row: the left two columns are for a model run on all plots, and the right two columns are for a model run on all 100% forest plots.Bottom row: the left two columns are for a model run on all plots with positive net biomass growth only, and the right two columns are for a model run on 100% forest plots with positive growth only.

Figure 8 .
Figure 8. Results of the best models for forest net biomass growth (NBG).Top row: the left two columns are for a model run on all plots, and the right two columns are for a model run on all 100% forest plots.Bottom row: the left two columns are for a model run on all plots with positive net biomass growth only, and the right two columns are for a model run on 100% forest plots with positive growth only.

Table 1 .
Image dates used to generate time series data and phenology products.

Table 3 .
Results for forest biomass models in the various plot groups.Adj.R 2 is the coefficient of determination adjusted for the number of predictors, carea means crown area, and Ht refers to the respective model's height variable.Phenology variables are as defined in Section 2.2.2.

Table 4 .
Results for forest net biomass growth models in the various plot groups.