Trends of land surface phenology derived from passive microwave and optical remote sensing systems and associated drivers across the dry tropics 1992–2012

Changes in vegetation phenology are among the most sensitive biological responses to global change. While land surface phenological changes in the Northern Hemisphere have been extensively studied from the widely used long-term AVHRR (Advanced Very High Resolution Radiometer) data, current knowledge on land surface phenological trends and the associated drivers remains uncertain for the tropics. This uncertainty is partly due to the well-known challenges of applying satellite-derived vegetation indices from the optical domain in areas prone to frequent cloud cover. The long-term vegetation optical depth (VOD) product from satellite passive microwaves features less sensitivity to atmospheric perturbations and measures different vegetation traits and functioning as compared to optical sensors. VOD thereby provides an independent and complementary data source for studying land surface phenology and here we performed a combined analysis of the VOD and AVHRR NDVI (Normalized Difference Vegetation Index) datasets for the dry tropics (25°N to 25°S) during 1992–2012. We find a general delay in the VOD derived start of season (SOS) and end of season (EOS) as compared to NDVI derived metrics, however with clear differences among land cover and continents. Pixels characterized by significant phenological trends ( P < 0.05) account for up to 20% of the study area for each phenological metric of NDVI and VOD, with large spatial difference between the two sensor systems. About 50% of the pixels studied show significant phenological changes in either VOD or NDVI metrics. Drivers of phenological changes were assessed for pixels of high agreement between VOD and NDVI phenological metrics (serving as a means of reducing noise-related uncertainty). We find rainfall variability and woody vegetation change to be the main forcing variables of phenological trends for most of the dry tropical biomes, while fire events and land cover change are recognized as second-order drivers. Taken together, our study provides new insights on land surface phenological changes and the associated drivers in the dry tropics, as based on the complementary long-term data sources of VOD and NDVI, sensitive to changes in vegetation water content and greenness, respectively.


Introduction
Seasonal variations in vegetation properties such as leaf area, photosynthetic productivity and water content are tightly coupled to global carbon, water, and energy fluxes (Bonan, 2008;Piao et al., 2007) and are recognized as a key indicator for ecosystem structure and functioning (Peñuelas et al., 2009).During recent decades, climate change has profoundly affected vegetation phenology across the globe through global warming and altered rainfall seasonality (Peñuelas and Filella, 2001) and altered phenophases in turn also affect climate via biophysical feedbacks (Richardson et al., 2013;Zeng et al., 2018).
Considerable efforts have been made on mapping phenology based on the combined use of ground observations and remotely sensed data (Linderholm, 2006), using local in situ measurements to train and validate remote sensing and modeling approaches at larger spatial scales (White et al., 2009;Xu et al., 2016).Periodical variations in vegetation as observed by remote sensing systems are usually termed land surface phenology in contrast to the term vegetation phenology that refers to life cycle events (e.g., bud break, flowering, and leaf senescence) of individual plants or species at the field level.Since collection of ground observations is laborious and limited in spatial extent, continuous satellite observations have become a widely recognized tool for monitoring large-scale vegetation phenology from regional to global scales.Based on remote sensing vegetation indices (VIs), such as Normalized Difference Vegetation Index (NDVI), land surface phenology is usually described by some key metrics i.e., start and end of growing season (SOS and EOS, respectively), from which phenological trends over time can be quantified (Verger et al., 2016;Zhang et al., 2003).
Long-term global land surface phenological trends have been extensively studied using satellite observations from the AVHRR (Advanced Very High Resolution Radiometer) sensors covering the period from 1981 to present.The AVHRR VI data measure the greenness of the canopy, which is closely linked to the leaf area, chlorophyll abundance and thus the overall vegetation productivity (Myneni and Hall, 1995).Particularly, for the Northern Hemisphere, several studies have shown widespread phenological changes in AVHRR VI based SOS and EOS (Garonna et al., 2014;Jeong et al., 2011;Piao et al., 2006), mainly as a consequence of global warming (Jeong et al., 2011).Contrastingly, studies of long-term phenological changes in the tropics are rather limited (Garonna et al., 2016).The few AVHRR-based studies conducted at a global scale revealed fewer areas of significant changes and more spatially heterogeneous land surface phenological change patterns in the tropics as compared to the Northern Hemisphere (Buitenwerf et al., 2015;Garonna et al., 2016).These differences in patterns of land surface phenology between tropics and the Northern Hemisphere are likely due to the high inter-annual variation in tropical land surface phenology related to climate oscillations (Brown et al., 2010) and complex relationships between driving factors including both climatic and anthropogenic (land management) forcing (Archibald and Scholes, 2007;Broich et al., 2014).Moreover, cloud-cover and atmospheric disturbances largely reduce the signal-to-noise level and availability of high-quality optical satellite imageries in cloud-prone tropical regions (Fensholt et al., 2007), introducing large uncertainties in the derived phenological metrics and the associated inter-annual changes (Ganguly et al., 2010;Zhang et al., 2009).
The vegetation optical depth (VOD) retrieved from satellite passive microwaves provides an independent and complementary data source to assess vegetation phenology at regional to global scales (Jones et al., 2011).VOD is linearly correlated with vegetation water content (Van de Griend and Wigneron, 2004;Jackson and Schmugge, 1991) and is sensitive to both leafy and woody vegetation components (Guglielmetti et al., 2007;Jones et al., 2013;Tian et al., 2017).VOD is only minimally affected by cloud and aerosol contaminations and passive microwave signals are also not impacted by variations in illumination conditions as being the case for VI's derived from optical/near infrared spectral bands.These characteristics of VOD have recently evoked substantial interests for global vegetation studies, including ecosystemscale plant hydraulics (Konings and Gentine, 2016;Tian et al., 2018), carbon stocks and dynamics (Brandt et al., 2018;Liu et al., 2015), woody/forest cover (Brandt et al., 2017a;van Marle et al., 2016), vegetation productivity (Teubner et al., 2018), crop yield estimation (Alemu and Henebry, 2017b), and land surface phenology (Jones et al., 2012).In particular, the long-term VOD dataset  provides an unprecedented opportunity for studying global changes in land surface phenology over a period comparable to the commonly used AVHRR VI datasets (Liu et al., 2011).The VOD product has data gaps under frozen soil conditions, thus hampering the calculation of phenological metrics such as SOS and EOS in large parts of the Northern Hemisphere, which is however not a problem in the tropics where valid VOD retrievals are available throughout the year.
It is thus of high relevance to perform a combined analysis of land surface phenological trends observed from the VOD time series and the AVHRR NDVI data records during 1992-2012 (the overlapping period of data availability) in tropical regions (25°N to 25°S) to gain more comprehensive insights into recent changes in the phenology of dry tropical ecosystems.In this study, we: (1) analyze agreement/differences in spatio-temporal patterns of VOD and NDVI derived phenological metrics of SOS and EOS for different land cover types and subcontinental regions and (2) examine potential drivers of VOD/NDVI observed phenological trends.

Data and methods
The overall data processing and analysis workflow is illustrated in Fig. 1 and detailed descriptions of the individual steps are provided below in Sections 2.1 to 2.6.

VOD data
Satellite passive microwave radiometers sense the thermal emissions from land surface quantified as brightness temperature (TB, the product of emissivity and physical temperature), from which VOD (dimensionless) can be retrieved.VOD is a measure of the strength of vegetation attenuation effect on the microwave signals, including both the leafy and woody components of the canopy.VOD was found to be linearly correlated with vegetation water content (kg m −2 ) and the relationship varies by microwave frequency (Wigneron et al., 2017).Lband passive microwaves (1-2 GHz) have a strong penetration capability (i.e., a relatively weak vegetation attenuation effect) due to the long wavelength (15-30 cm), while C-/X-/K-band passive microwaves (6-18 GHz) are weaker in penetrating through the vegetation layer (Santi et al., 2009).Therefore, for dense deciduous forests, L-band VOD mainly senses water stored in trucks and branches with a relatively narrower seasonal dynamic range; whereas, C-/X-/K-band VODs mainly represent the water content fluctuations at canopy/crown level (leaves and branches) with larger seasonal dynamics.VODs from all these bands are minimally affected by clouds and aerosols.However, open surface water bodies have a strong influence on the accuracy of retrieved VOD, which needs to be carefully dealt with by masking out pixels accordingly.
The long-term VOD dataset employed here was produced and provided by Liu et al. (2011Liu et al. ( , 2015) ) at a spatial resolution of 0.25°and with a daily temporal resolution for the period 1992-2012.The VOD dataset was retrieved from several satellite microwave radiometers at frequencies ranging from 6.8 GHz to 19.4 GHz (C-/X-/K-band), including the Special Sensor Microwave Imager (SSM/I), the Advanced Microwave Scanning Radiometer for EOS (AMSR-E), the WindSat, and the FengYun-3B, using the land parameter retrieval model (Owe et al., 2001).A cumulative distribution function matching approach was applied to merge the retrievals from different sensors to form a consistent long-term VOD time series, having a dynamic range from 0 to 1.2 (Liu et al., 2011).For areas dominated by herbaceous vegetation, the seasonal variation in the VOD signal reflects the annual cycle of the entire herbaceous layer.For areas dominated by woody vegetation, the VOD seasonality is related to the canopy water status, i.e., fluctuations in the canopy leaf area and the associated water content in both leaves and branches (Momen et al., 2017;Santi et al., 2009;Tian et al., 2016).As the VOD seasonality may be influenced by seasonal inundation, we masked out pixels with a spatial extent of inundation > 10% (Supplementary Fig. S1a), based on a global inundation map (Fluet-Chouinard et al., 2015).

AVHRR NDVI data
We used the Global Inventory Modeling and Mapping Studies (GIMMS) NDVI 3rd generation version 1 product (Pinzon and Tucker, 2014) due to its high temporal consistency across multiple AVHRR sensors (Tian et al., 2015) and its frequent use for long-term phenological trend analysis.This NDVI dataset is provided in a 1/12°spatial and a semi-monthly temporal resolution by compositing (selecting the maximum value) the daily observations within each 15-day window to reduce the impacts from cloud cover and atmospheric influence.We aggregated the NDVI data to a 0.25°by pixel averaging to match the spatial resolution of the VOD dataset.

Extraction of phenological metrics
We applied the same method for VOD and NDVI time series to extract SOS and EOS, although VOD and NDVI respond to different biophysical characteristics (vegetation water content vs greenness) and show different seasonal patterns depending on ecosystem functional types (Fig. 2).SOS and EOS extractions were done at the per-pixel level using the TIMESAT software (Jönsson and Eklundh, 2004) by calculating the day of year (DOY) when VOD/NDVI reach half of their seasonal amplitudes (the difference between annual minimum and maximum) before and after the peaking time, respectively (Garonna et al., 2014;White et al., 2009).These settings are comparable to the middle point method applied by Garonna et al. (2016) used for an NDVI-based global phenological trend analysis.
The VOD and NDVI raw time series were smoothed using a Savitzky-Golay filter (window size = 60 days, i.e., 61 observations for VOD and 5 for NDVI) in the TIMESAT software (Fig. 2).Areas characterized by a seasonal amplitude < 0.1 (Supplementary Fig. S1b) and with an annual mean value < 0.1 (Supplementary Fig. S1c) for either VOD or NDVI were masked out, following the criteria used by Garonna et al. (2016).We also masked out areas with more than one growing season (Supplementary Fig. S1d) based on the phenology product generated by the Vegetation Index and Phenology Lab (http://www.vip.arizona.edu).Rather than reducing the VOD temporal resolution to match with NDVI (15-day), we used the daily VOD as input to obtain the highest possible accuracy of SOS and EOS retrievals to support our study aims.An intercomparison of SOS and EOS estimates as based on time series of daily VOD and 15-day composited VOD was conducted (Supplementary Fig. S2) showing high consistency (R 2 values and RMSE) between output results.

Comparison between VOD and NDVI phenological metrics
We examined the per-pixel intra-annual difference of each phenological metric (SOS and EOS) between VOD and NDVI.We also compared the differences in VOD and NDVI long-term phenological trends (SOS and EOS), using the non-parametric Theil-Sen linear trend analysis and the Mann-Kendall significance test (Mann, 1945;Sen, 1968).These intra-and inter-annual phenological differences were compared based on land cover types (using the European Space Agency (ESA) Climate Change Initiative (CCI) products, https://www.esa-landcovercci.org/).We regrouped the original land cover classes into grassland, shrubland, cropland, mosaic tree/shrub/grass/crop, and broadleaf deciduous forest (Fig. 3; pixels classified as broadleaf evergreen forest were excluded due to the limited seasonality producing high uncertainty in the SOS and EOS extraction), each representing different characteristics of ecosystem structure and functioning.Yet, vegetation characterized by the same land cover type may show a different biophysical behavior across space; for example, the rain use efficiency of the African semi-arid regions varies significantly between the Sahel and southern Africa (Ratzmann et al., 2016).Therefore, we also broke down the analysis to a sub-continental level by studying four individual regions, that is, northern sub-Saharan Africa (N.Africa), southern Africa (S.Africa), eastern South America (S.America), and northern Australia (N.Australia) (Fig. 3).
We further evaluated whether two independent groups were equivalent or not in the SOS and EOS difference as derived from VOD and NDVI.Test of significant difference between groups of variables measured across space from remote sensing can be ambiguous due to a large number of pixels and spatial autocorrelation of these pixels (Foody, 2009a(Foody, , 2009b)).To avoid this issue, we instead tested for equivalence between groups as has been done by Robinson andFroese (2004) andde Beurs et al. (2015).Specifically, we tested at a certain level (α) if the mean value of two groups were equivalent or not, with the null hypothesis: By rejecting H 0 , we conclude that the two tested groups are equivalent (the difference between them is not sufficient to care about).We set α to 25% of the mean standard deviation (Robinson and Froese, 2004) of the observed differences for all the tested groups (different land cover types and geographic regions) corresponding to 5.26 days.

Assessment of NDVI and VOD data quality by inter-comparison
Both NDVI and VOD time series data are subject to some level of noise and perturbation effects, especially for the humid and densely vegetated areas.NDVI suffers from atmospheric contamination from aerosols, water vapor and extensive cloud cover (Vermote et al., 2002) and directional effects from varying sun-sensor geometry (Fensholt et al., 2010;Leeuwen and Orr, 2006), whereas VOD retrieval accuracy decreases when vegetation exceeds a certain level (Jones et al., 2011).Also, both NDVI and X-band VOD show saturation effects as a function of vegetation density (Jones et al., 2011;Sellers, 1985;Tian et al., 2016).These well-known issues are likely to introduce uncertainties in the calculated phenological metrics and will thereby impact the interannual variations of the respective metrics.VOD and NDVI are expected to co-vary between successive years in response to climate factors (e.g. earlier rains in a specific year will result in earlier SOS for both VOD and NDVI), except for years with abrupt disturbances.However, this short-term co-variation will not necessarily result in identical long-term trends if changes in structure and/or species composition are occurring.We computed the Pearson's product-moment correlation coefficient between detrended VOD and NDVI phenological metrics (SOS and EOS) during 1992-2012 to assess the agreement between these two independent datasets.Detrended time series of phenological metrics were used to focus this part of the analysis on the short-term co-variation between VOD and NDVI.Pixels with significant (P < 0.05) correlation for either SOS or EOS were marked as "high agreement" between VOD and NDVI (good data quality), otherwise marked as "low agreement".
As VOD and NDVI measure different vegetation traits (water content and greenness, respectively) it should be noted that a per-pixel low correlation coefficient could potentially also be caused by changes in species composition (e.g., woody vegetation fraction), climatic variation (e.g., onset of rainy season) and land cover disturbances (e.g., fires) because the response of VOD and NDVI to the above-mentioned changes might not be identical.However, the aim here was to make a first assessment of the importance of impacts from various perturbations on both VOD and NDVI derived phenological metrics, to be used in our following analysis of exploring potential drivers of phenological trends.

Examining potential drivers of phenological trends
A significant trend in phenology is usually associated with structural changes in vegetation composition and/or functioning (Peñuelas et al., 2009).Here, we label a pixel to be characterized by a change in phenology if trends in any of the four examined metrics (VOD SOS, VOD EOS, NDVI SOS, and NDVI EOS) were significant.We subsequently examined the potential drivers of a significant phenological change by analyzing various driver variables.For a given land cover type at the sub-continental level, we compared the group of pixels showing significant phenological trends against the reference group of pixels with non-significant phenological trends (Fig. 1).A given driver variable was identified as potential driver if its value was not equivalent between the groups of significant and non-significant phenological trends (see Section 2.4).A two-sided equivalent test was used for all the comparisons except for the analysis of land cover change intensity as a potential driver, where a one-sided test was used with the null hypothesis: We selected four potential controlling variables, which were land cover change and temporal dynamics of rainfall, fire events, and woody vegetation (encroachment or deforestation/forest degradation).Here, temporal dynamics refer to the per-pixel long-term trend (slope of linear regression against time) and inter-annual variation (IAV, the anomaly to the fitted trend line).These drivers may not be mutually exclusive, but rather represent complementary information about forcing mechanisms leading to a change in vegetation composition and/or functioning.For example, a shrubland area characterized by woody encroachment may show altered ecosystem functioning even though the land cover remains categorized as shrubland, but eventually the   S1a-d.land cover might change as well (e.g.into savanna woodland).Only pixels assigned with high agreement between VOD and NDVI phenological metrics (see Section 2.5) were included in this analysis to reduce the potential impacts from data artifacts (atmospheric perturbations, viewing geometry, signal saturation etc.) on the attribution of drivers.
The prevalence of land cover change (the number of land cover types of a given location during the period of analysis, here referred to as land cover change intensity) was calculated based on the annual ESA CCI land cover time series.We first counted the number of land cover types during 1992-2012 for each pixel at the original spatial resolution (300 m) and then resampled the output to VOD pixel size by averaging.Annual rainfall data from CHIRPS v2.0 (Funk et al., 2015) covering the period 1992-2012 were used and the data provided at a spatial resolution of 0.05°were resampled to 0.25°by averaging.Changes in woody vegetation (trees and shrubs) estimated using annual minimum VOD observations (corresponding to dry season observations with no/limited herbaceous cover), found to be a robust proxy for fractional woody cover (Brandt et al., 2017a).This VOD dataset is found to be homogeneous across sensors as assessed by Tian et al. (2016).Data on fire events were obtained from the GFED4 monthly burned area product (excluding small fires) during 1996-2012 (four years less than the other datasets) at 0.25°spatial resolution (Giglio et al., 2013).Per-pixel monthly burned area information was aggregated to annual sums for further analysis.

Differences between VOD and NDVI phenological metrics
Mean values for SOS and EOS of VOD generally occurred later in the season than those derived from the NDVI time series  for the majority of pixels analyzed (83% and 81% of SOS and EOS for the study area, respectively) (Fig. 4 and Tables 1-2, see Supplementary Fig. S3 for the mean DOY and Supplementary Fig. S4 for the standard deviation of VOD and NDVI phenology difference).A larger difference between VOD and NDVI is seen for SOS as compared to EOS in most areas and VOD and NDVI phenological differences also show clear spatial patterns clustered within different sub-continental regions.In the N. Africa subset, the southernmost Sudan-Guinean zone is characterized by a higher vegetation density (deciduous forests and shrubs) with larger differences in SOS (Fig. 4a and Table 1) and EOS (Fig. 4b and Table 2) than the Sahelian zone to the north (dominated by grassland and shrubland).A similar cluster of differences between VOD and NDVI co-varying with vegetation density was also found in southern Africa for SOS, but being less clear for EOS.The analysis revealed that the same land cover type exhibits significantly different patterns in different subcontinental regions (Tables 1-2).For example, a widespread earlier SOS of VOD as compared to NDVI was found in the N. Australian shrublands, which was not found to be the case for shrublands in the other tropical regions.

Inter-annual trends in VOD and NDVI phenological metrics
Overall, the percentage of pixels characterized by significant trends (P < 0.05) is covering up to 20% of the study area for each category of SOS and EOS as derived from both VOD and NDVI, with spatial clusters showing both significantly positive and negative trends (Fig. 5).Significant trends towards later dates of both SOS and EOS for N. Africa region were observed from NDVI, which was less prominent in the analysis based on VOD data.For southern Africa, both VOD phenological metrics (SOS and EOS) show large spatially clustered patterns of significant trends towards both earlier and later dates, whereas the NDVI metrics show a more scattered pattern.In S. America notably, the overall EOS trend in the Brazilian Cerrado area shows opposite patterns between VOD (later) and NDVI (earlier).
When combining the results of VOD and NDVI trends, both SOS and EOS show significant trends for around 25% of the study area, each with the majority being characterized by only one significant metric for either VOD or NDVI, while few pixels show converging or diverging VOD/NDVI significant trends for SOS and EOS (Fig. 6a, b).Detailed joint occurrence tables showing all combinations of VOD and NDVI trends (i.e., significant positive, significant negative and non-significant) for each sub-continental region and land cover type are provided as Supplementary material (Tables S1-S18).Generally, the N. Africa region show more significant phenological trends in NDVI, while other regions show more or less equivalent proportions of significant trends between VOD and NDVI.A significant trend in any of the four metrics examined (VOD SOS, VOD EOS, NDVI SOS, and NDVI EOS) exists for about 50% of the area studied (Fig. 6c).

Agreement between VOD and NDVI phenological metrics
The correlation between inter-annual variations in VOD and NDVI metrics (both SOS and EOS; Fig. 7a, b) are higher in the drier regions, e.g., the African Sahel, grasslands of southern Africa, Brazilian Caatinga, and northern Australia as compared to the wetter regions, e.g., the African Sudan-Guinean, the wet Miombo woodland, north-eastern Cerrado Brazil, and southeast Asia.These drier regions correspond well with previously reported areas where NDVI data have ample high quality cloud-free observations (Fensholt and Proud, 2012).When combining results from both SOS and EOS (Fig. 7c), about 70% of the study area shows high agreement (significant correlation, P < 0.05) between VOD and NDVI phenological metrics (from where the following analysis on potential drivers was performed).It is evident that trend agreement is lower for the S. America continental subset as compared to the three other subsets.

Potential drivers of observed phenological trends
Table 3 shows results of the equivalence tests between pixels with and without significant phenological trends (Fig. 6c) for each sub-continental region and land cover type, while the drivers identified are summarized in Table 4. Overall, the temporal dynamics in woody vegetation, rainfall, and fire events were all important drivers for the study area whereas land cover change was found to be less dominant.Whereas the variables examined had little explanatory power as drivers of phenological trends for land cover types of N. Africa, all drivers collectively were found to play a role for phenological changes of land cover types in S. America.For southern Africa, woody vegetation dynamics play a major role for all the land cover types (except for the deciduous forest areas), coinciding with areas documented by extensive woody encroachment (O'Connor et al., 2014).For N. Australia, the analysis revealed dynamics in woody vegetation and fire events as drivers of phenological trends for the land cover types of grassland and mosaic tree/shrub/grass/crop.Cropland and broadleaf deciduous forest were excluded from the sub-continental analysis of N. Australia due to the limited number of pixels available (< 50).The summary of identified drivers across land cover types and continents (Table 4) show that rainfall and woody vegetation dynamics are the most frequently occurring variables explaining phenological trends (nine occurrences each), followed by fire events (seven occurrences) and land cover change intensity (two occurrences).

Biophysical comparability of VOD and NDVI phenology
Similar to NDVI, the VOD product employed in this study also shows a clear seasonal response to vegetation growth (Fig. 2), which allows for VOD-based mapping of land surface phenology in the same way as NDVI (Jones et al., 2011;Jones et al., 2012).The observed land surface phenology in both NDVI and VOD is driven by vegetation seasonal variations but monitoring of phenology from optical and microwave sensor systems involves different sensitivity to biophysical traits of vegetation, which complicates the interpretation of observed patterns in phenological metrics and changes herein.While NDVI measures the vegetation greenness related to the leaf component, VOD is sensitive to the vegetation water content in both the leaf and woody components.Therefore, the phenological difference derived from NDVI and VOD should be interpreted as the sum of both the leaf water status at a given level of greenness and the water status of the woody vegetation components.At a seasonal scale, canopy greenness and leaf water content are usually highly correlated, both co-varying with the amount of green vegetation biomass.This covariation is supported by the small SOS/EOS difference between VOD and NDVI for the grasslands (Figs. 2 and 3) and is also consistent with a previous study focusing on homogenous croplands of northern Eurasia (Alemu and Henebry, 2017a).In contrast, the seasonal development of water content of the woody component does not necessarily follow canopy greenness, often exhibiting a certain time lag (Borchert, 1994;Tian et al., 2018), which explains the larger SOS/EOS difference for the forested areas (Figs. 2 and 3).For woodlands, the NDVI increases rapidly and reaches the highest level (around 0.8) following the green-up period and then remains at this level (or even decreases a bit) during the growing season.In the event of NDVI saturation effects, the estimated SOS will be underestimated (detected earlier), but the underestimation should be far below the observed SOS difference between VOD and NDVI (around 2 months in woodlands).Overall, the observed delay in SOS and EOS from VOD as compared to NDVI indicate general asynchronicity in the vegetation seasonal cycle of canopy greenness and vegetation water content.However, the differences between VOD and NDVI metrics show sub-continental region-specific patterns for the same land cover type (Fig. 4), suggesting the existence of distinct eco-physiological traits of vegetation with regionspecific functional and structural adaptation to biotic and abiotic factors.
Since the VOD product used in this study is retrieved from relatively high microwave frequencies (> 6 GHz), the variation in vegetation water content mainly origins from the canopy/crown level for dense forest areas.This is different for the VOD retrieved from SMOS (Soil Moisture and Ocean Salinity) at L-band (1.4 GHz), which can sense the water storage in stems of dense forest (Brandt et al., 2018;Tian et al., 2018).However, SMOS VOD data are only available since 2010, and the short period and noisy signal limit its current suitability for assessing trends in land surface phenology.S1a-d.

Table 1
The mean number of days of difference between VOD SOS and NDVI SOS for each land cover type and sub-continental region.The capital letters (e.g., A, B) indicate the equivalence test between geographic regions for a given land cover type (individual columns), while the numbers (e.g., 1, 2) indicate the equivalence test between land cover types for a given geographic region (individual rows).Groups with the same capital letter or number are equivalent in values.The equivalence level α is 5.26 days.tropical areas (Fensholt et al., 2007;Rankine et al., 2017).Temporal compositing (e.g., semi-monthly or monthly maximums) or smoothing (gap-filling) is often used to reduce the impact from atmospheric contamination in NDVI time series (Holben, 1986;Zhang, 2015).The 15day temporal resolution of the GIMMS AVHRR NDVI inevitably will introduce uncertainties in the extracted phenological metrics (i.e., SOS and EOS) as compared to analyzes conducted on daily data.Also, the coarse spatial resolution may conceal details in the spatial domain (Zhang et al., 2017) with uneven levels of uncertainty between SOS and EOS (Zhang et al., 2018).Yet, the GIMMS AVHRR NDVI dataset has been widely used in detecting phenological changes, particularly in the Northern Hemisphere (Garonna et al., 2014;Jeong et al., 2011;Piao et al., 2006;Xu et al., 2016).The changes detected from AVHRR time series spanning several decades are commonly reported to be at a rate of a few days, which is a change rate being much shorter than the temporal resolution of the composited dataset, as summarized by Garonna et al. (2016).The results are however valid, assuming that data uncertainties are randomly distributed over long periods, and patterns of change are prevalent for the Northern Hemisphere showing noticeable temperature-driven seasonal trends (Zhang et al., 2009;Zhang et al., 2004).Contrastingly, phenological changes in the dry tropics are often more  Light grey colored areas are masked according to information provided in Supplementary Fig. S1a-d.
spatially heterogeneous as potentially driven by more than one dominant variable (Garonna et al., 2016).Our study provides the first spatially comprehensive evaluation of GIMMS AVHRR NDVI phenological trends in the dry tropics by comparing with independent results from VOD data.It is by nature a difficult task to assess how severe is the impact from NDVI and VOD data uncertainties on phenological trend retrievals, as differences between metrics from NDVI and VOD can be caused by both data quality and the different sensitivity to greenness and canopy water content, combined with impacts from possible changes in climate and land management.However, areas of low correlation between detrended NDVI and VOD phenological metrics matches well with areas of widespread cloud cover, thus suggesting a cautious interpretation of the NDVI phenological trends reported here (Fig. 5) and in previous studies (Buitenwerf et al., 2015;Garonna et al., 2016).We attributed the potential drivers only within areas of the dry tropics showing high SOS/EOS agreement between NDVI and VOD, which reduces the effects of data artifacts and significant uncertainties.Hence, we did not further explore the phenological trends and potential drivers from VOD data in areas of no correlation with NDVI and in areas masked out based on a common NDVI/VOD seasonal amplitude threshold (here masked based on a threshold of 0.1).However, as VOD is less affected by atmospheric contamination (even cloud cover) and shows less saturation effects as compared to NDVI (Tian et al., 2016), VOD-based information is likely to contain interesting information on land surface phenology in the areas excluded here.Future studies on land surface phenology in these cloud-prone tropical regions should also explore the potential use of data from geostationary satellite systems (e.g. the Spinning Enhanced Visible and Infrared Imager carried by Meteosat satellites) producing more cloud-free observations as compared the polarorbiting sensors with implications for phenology studies in the humid rainforests (Fensholt et al., 2007;Fensholt et al., 2006;Guan et al., 2014;Julien et al., 2012;Yan et al., 2016a;Yan et al., 2017;Yan et al., 2016b).

Dry tropical phenological trends
The land surface phenological trends of SOS assessed by NDVI during 1992-2012 (this study) show visually different spatial patterns with those reported by Garonna et al. (2016), who used the same dataset and methods but during 1982-2012.With a 10-year shorter period, our study shows visually more areas characterized by significant trends in SOS, for example, the significant delay in the N. African Fig. 6.Comparison of (a) SOS trends and (b) EOS trends between VOD and NDVI and (c) the combined mapping of significant trends in any of the phenological metrics.Pixels with significant trends (P < 0.05) are colored, while non-significant trends are shown in dark grey.The per-category areal percentages for each selected sub-continental region (N: N. Africa, S: southern Africa, SA: S. America, A: northern Australia) and land cover type (G: grassland, S: shrubland, C: cropland, M: mosaic tree/shrub/grass/crop, D: broadleaf deciduous forest) are inserted as bar-plots.The width of bars indicates the relative size of each region and land cover type (wider means a larger areal size).Light grey colored areas are masked according to information provided in Supplementary Fig. S1a-d.drylands (Fig. 5b).Yet, also similar EOS trend patterns were observed between these two studies, for instance, the widespread significant delaying trend in N. Africa and the significant earlier trend in Brazil (Fig. 5d).The difference in trends obtained between different periods is likely related to the vegetation response to climate variations/perturbations as herbaceous and woody vegetation in dry tropics are sensitive to water availability (Huber et al., 2011) and some species of woody vegetation can quickly die off due to a lack of rainfall in several successive years (Brandt et al., 2017b).Contrastingly, similar trend patterns during 1992-2012 and 1982-2012 in some areas indicate relatively persistent shifts in vegetation structure/species composition and/ or climate regimes.
The large discrepancy between VOD and NDVI phenological trend patterns in both SOS and EOS (Fig. 6a, b) shows that ecosystem-scale vegetation changes may not always translate into uni-directional changes in vegetation functional traits, highlighting the benefits of including independent and complementary data sources for studying land surface phenology.The combined use of VOD and NDVI can thus be used for a more comprehensively understanding of phenological changes as compared to using only a single dataset.For example, the Sahelian grasslands show similar phenology dates between VOD and NDVI dataset (Fig. 4) but inconsistent long-term trends (Fig. 7), which could be an indication of changes in herbaceous species (e.g., a change from horizontal to vertical leaf orientation may induce larger changes in NDVI (Mbow et al., 2013) than in VOD).Furthermore, VOD and NDVI even show some scattered areas of significant diverging phenological trends (notably in southern Africa and S. America; Fig. 6a,  b and Supplementary Tables S12-S13), which could be caused by extensive land transformations such as large-scale deforestation caused by cropland expansion in the Brazilian Cerrado (Lapola et al., 2013).

Drivers and implications of observed phenological trends
Understanding drivers of observed trends in land surface phenology across the dry tropics is important but also challenging.Rather than analyzing the per-pixel inter-annual development of land surface phenological metrics and potential drivers, (i.e., land cover change and dynamics in rainfall, fire, and woody vegetation), we here opted for a spatial approach based on spatial statistics of the potential drivers within the same land cover type and geographic region (Tables 3-4).Our results provide evidence for the effects of trends and variations (IAV) in woody vegetation, rainfall, and fire events on phenological changes over the dry tropics.Woody encroachment was reported to be widespread across the tropical drylands (Brandt et al., 2016;Brandt et al., 2018;O'Connor et al., 2014;Stevens et al., 2017) and significant negative trends in fire frequency have been documented for the tropics (Andela et al., 2017;Andela and van der Werf, 2014), which are both consistent with the widespread phenological changes found in our study.
Recent studies documented the widespread pre-rain green-up phenomenon in African deciduous forests (Ryan et al., 2017;Tracy et al., 2018), which suggest a decoupling of SOS of woody vegetation and the Light grey colored areas are masked according to information provided in Supplementary Fig. S1a-d.onset of the rainy season.The controlling process of the pre-rain greenup is still unclear (Stock, 2017), although different environmental cues have been suggested to trigger the leaf emergence of woody vegetation in these areas, including air humidity (Alemu and Henebry, 2017b;Brown and de Beurs, 2008;Do et al., 2005) and photoperiod (Ryan et al., 2017).In any case, the vegetation water consumption before rainfall onset is expected to come from ground water which is recharged from previous years' rainfall.Therefore, these facts do not contradict the dominant role of rainfall variability on phenological trends for the southern African broadleaf deciduous forests as revealed by this study (Tables 3-4).Moreover, the magnitude of leaf development during the pre-rain green-up period is usually smaller as compared to that during the rainy period (Ryan et al., 2017).Also, the specific phenological metrics calculated here were based on the middle point between the annual minimum and maximum VOD/NDVI values, which will consequently reduce the effect of pre-rain green-up from the trees while enhance the sensitivity to the rainfall-driven herbaceous layer.Our results highlight the importance of rainfall variability in the tropical land surface phenology, supporting the hypothesis that vegetation structure of seasonally dry tropical forests will be sensitive to present and future changes in rainfall regimes (Kara et al., 2017).Therefore, it is likely that a projected shift in rainfall regimes due to climate change (Feng et al., 2013;Seth et al., 2013) may introduce further changes in land surface phenology, as a consequence of changes in the coexistence between herbaceous and woody vegetation in dry tropical forest ecosystems (Zhang et al., 2019).
Alterations in ecosystem functioning caused by land cover change are expected to impact vegetation in the form of an abrupt shift, after which the ecosystems would be still predominantly controlled by other factors, such as rainfall, in the tropical drylands.This mechanism could explain the limited impacts of land cover change in driving land surface phenological trends.Areas without any identified drivers for N. African (grassland, shrubland, and broadleaf deciduous forest) and N. Australia (shrubland, Tables 3-4) indicate other forcing variables not included in the current analysis (e.g.grazing, temperature and radiation).

Table 3
Mean values of examined drivers for pixels characterized by non-significant (NS) and significant (S) phenological trends, respectively, for each land cover type and sub-continental region.Presence of a phenological trend was defined by considering both VOD and NDVI metrics (Fig. 6c) and only pixels with high agreement (Fig. 7c) were analyzed.The equivalence test was performed between the groups of NS and S for each land cover type and sub-continental region.Numbers in bold within each of the driver variables indicate that groups are not equivalent (i.e. the difference between them is significant) at the respective level of α and thereby regarded as drivers of significant phenological trends.

Table 4
Summary of identified drivers of significant phenological trends for each land cover type and sub-continental region.The results were combined for the equivalence tests of trend and interannual variability (IAV) of woody vegetation, rainfall, fire events (burnt area) and land cover change intensity (LCC), respectively.

Conclusion
We analyzed land surface phenological trends in the dry tropical regions (25°N to 25°S) during 1992-2012 by a combined analysis of VOD and AVHRR NDVI datasets based on passive microwave and optical remote sensing systems, respectively.The findings represent a new way of mapping and understanding recent changes in phenology of the dry tropical ecosystems building on the complementarity of the data sources used, being sensitive to different vegetation functional traits and characterized by different abilities to overcome noise-propagation in the phenology retrievals from atmospheric perturbations and saturation effects.We found a general delay of the VOD derived start of season (SOS) and end of season (EOS) as compared to NDVI derived metrics.Pixels characterized by significant phenological trends covered up to 20% of the study area for each phenological metric (SOS/EOS from VOD/NDVI), however showing large spatial differences in areas of significant trends between the two sensor systems.Temporal dynamics in rainfall, woody vegetation, followed by fire events and land cover change were identified to be the drivers of phenological trends (in decreasing order of importance), however varying across land cover types and sub-continental regions.Further research should be specifically designed to take advantage of the complementary biophysical bearing of VOD (sensitive to vegetation water content) and NDVI (sensitive to canopy greenness) for studies of land surface phenological changes and associated ecophysiology in the tropics, where the single use of optical sensor systems is often hampered by prevailing cloud cover.

Fig. 1 .
Fig. 1.Flow chart of the study involving extraction of VOD and NDVI phenological metrics, inter-comparison and analysis of drivers of phenological changes.

Fig. 2 .
Fig. 2. Examples of VOD and NDVI seasonal variations during 2005 at (a) a grassland site (location: 15.25N, 15W) in northern Senegal and (b) a woodland site (location: 6.75N, 17.25E) in the southern Central African Republic.Crosses represent the VOD and NDVI raw values and lines represent smoothed seasonal curves using the Savitzky-Golay filter.The solid points represent retrieved start and end of growing season as based on VOD and NDVI.The daily rainfall data (dark grey bars) are from the CHIRPS rainfall product.

Fig. 3 .
Fig. 3. Land cover types of the dry tropics (25°N to 25°S), based on the ESA CCI land cover class data in 2000 (representing the middle of our study period).The four subcontinental regions examined are indicated by blank bounding boxes.Light grey colored areas are masked according to information provided in Supplementary Fig. S1a-d.
NDVI phenologyNDVI is influenced by atmospheric contamination, and notably cloud cover largely reduces the number of valid NDVI observations in

Fig. 4 .
Fig. 4. Mean number of of difference between (a) VOD SOS and NDVI SOS and between (b) VOD EOS and NDVI EOS during 1992-2012.Light grey colored areas are masked according to information provided in Supplementary Fig. S1a-d.

Fig. 5 .
Fig. 5. Linear trends in phenological metrics ofVOD and NDVI (1992-2012)  for SOS (a, b) and EOS (c, d).Pixels with significant trends (P < 0.05) are colored, while non-significant trends are shown in dark grey.The areal percentages of the three types of trends (i.e., earlier, later, and non-significant) for each selected subcontinental region (N: N. Africa, S: southern Africa, SA: S. America, A: northern Australia) and land cover type (G: grassland, S: shrubland, C: cropland, M: mosaic tree/shrub/grass/crop, D: broadleaf deciduous forest) are inserted as bar-plots.The width of bars indicates the relative size of each region and land cover type (wider means a larger areal size).Light grey colored areas are masked according to information provided in Supplementary Fig.S1a-d.

Fig. 7 .
Fig. 7. Correlation between VOD and NDVI derived SOS (a) and EOS (b) as well as VOD and NDVI phenological metric agreement during 1992-2012 as evaluated by the existence of a significant correlation (P < 0.05) in either SOS or EOS (c).The areal percentages of each category for each selected sub-continental region (N: N. Africa, S: southern Africa, SA: S. America, A: northern Australia) and land cover type (G: grassland, S: shrubland, C: cropland, M: mosaic tree/shrub/grass/crop, D: broadleaf deciduous forest) are inserted as bar-plots.The width of bars indicates the relative size of each region and land cover type (wider means a larger areal size).Light grey colored areas are masked according to information provided in Supplementary Fig.S1a-d.

Table 2
The mean number of days of difference between VOD EOS and NDVI EOS for each land cover type and sub-continental region.The capital letters (e.g., A, B) indicate the equivalence test between geographic regions for a given land cover type (individual columns), while the numbers (e.g., 1, 2) indicate the equivalence test between land cover types for a given geographic region (individual rows).Groups with the same capital letter or number are equivalent in values.The equivalence level α is 5.26 days.