Response of Land Surface Phenology to Variation in Tree Cover during Green-Up and Senescence Periods in the Semi-Arid Savanna of Southern Africa

: Understanding the spatio-temporal dynamics of land surface phenology is important to understanding changes in landscape ecological processes of semi-arid savannas in Southern Africa. The aim of the study was to determine the inﬂuence of variation in tree cover percentage on land surface phenological response in the semi-arid savanna of Southern Africa. Various land surface phenological metrics for the green-up and senescing periods of the vegetation were retrieved from leaf index area (LAI) seasonal time series (2001 to 2015) maps for a study region in South Africa. Tree cover (%) data for 100 randomly selected polygons grouped into three tree cover classes, low (<20%, n = 44), medium (20–40%, n = 22) and high (>40%, n = 34), were used to determine the inﬂuence of varying tree cover (%) on the phenological metrics by means of the t -test. The differences in the means between tree cover classes were statistically signiﬁcant ( t -test p < 0.05) for the senescence period metrics but not for the green-up period metrics. The categorical data results were supported by regression results involving tree cover and the various phenological metrics, where tree cover (%) explained 40% of the variance in day of the year at end of growing season compared to 3% for the start of the growing season. An analysis of the impact of rainfall on the land surface phenological metrics showed that rainfall inﬂuences the green-up period metrics but not the senescence period metrics. Quantifying the contribution of tree cover to the day of the year at end of growing season could be important in the assessment of the spatial variability of a savanna ecological process such as the risk of ﬁre spread with time.


Introduction
Spaceborne remote sensing has provided a cost effective means to study, understand and monitor spatio-temporal changes in vegetated land surfaces [1].The study of the seasonal pattern of variation in vegetated land surface using remote sensing is known as land surface phenology [1,2], while the traditional phenology is the study of the recurring patterns of vegetation growth and development and their relation to climate, usually for single vegetation types or species [3][4][5].For example, the phenology of deciduous trees in the Southern African savanna is characterised by periods of leaf green-up, leaf display and leaf fall [6].On the other hand, the land surface phenological signature of a remotely sensed pixel in the African savanna is determined by the spectral characteristics of the various vegetation types present in the pixel, which could be any combination of types including deciduous trees, evergreen trees, perennial or annual grasses [1].
Commonly defined land surface phenological metrics include the dates corresponding to the start, end and length of the growing season.These metrics are important parameters in ecosystem simulation models and coupled biosphere/atmosphere models [7].For example, the timing and the length of the growing season control spatio-temporal variations of carbon and water cycles and influences latent/sensible heat transport [7,8].Wessels et al. [9] used land surface phenological metrics to describe and map various vegetation biomes in South Africa.Archibald and Scholes [6] have stressed the need for reliable phenological models with the advent global change given the influence of climatic factors on phenological metrics.The land surface phenology of Southern African semi-arid savannas is still largely understudied.
The savanna biome is characterised by a mixture of grasses (herbaceous vegetation in general) and trees of varying ratios e.g., more trees in general means less grass [10,11].The grass layer is made up of annual and perennial grass species [12,13], while the trees consist of deciduous and evergreen tree species of various heights and crown dimensions [14][15][16].Savanna grasses and trees have been shown to exhibit differential leaf phenological patterns because of differences in the factors that control the phenological events [6,17].In general, Southern African savanna phenology is influenced by climate, topography, grazing or browsing activities and fire [18,19].Grass production is generally limited by rainfall or water availability [6,10].The start of the growing season of the grass layer has been shown to be closely associated to rainfall [20].On the other hand, savanna tree species exhibit varying phenological patterns.For example, Whitecross et al. [20] recently demonstrated that the green-up dates amongst savanna tree species are more variable when compared to the differences between grass species.This finding supports the earlier study of Archibald and Scholes [6] who showed that most trees in the Southern part of the Kruger National Park, South Africa flushed new leaves before the onset of the rainy season with the exception Combretum apiculatum and C. Zeyheri which started greening up at the same time with grasses.With respect to interannual variability of the phenological metrics, Archibald and Scholes [6] demonstrated savanna trees show lower interannual variability in phenology in comparison to grasses.
Studying the phenology of savanna trees or grasses using commonly available satellite imageries such the Advanced Very High Resolution Radiometer (AVHRR) and the Moderate Resolution Imaging Spectroradiometer (MODIS) is not feasible because of the difficulties to separate trees and grasses using these coarse resolution sensors.AVHRR (1km resolution) or MODIS (250-500 m) pixels are most often mixed pixels of varying proportion of trees, grasses and bare soil [1].As described in Archibald and Scholes [6], various studies have made different assumptions in studying the phenology of trees and grasses using the above coarse resolution sensors and reached varying conclusions.For example, over-prediction of early green-up was observed in sites with high grass cover when the grass layer was not considered [21].In the domain of land surface phenology, the objective is not to study the phenology of trees and grass separately but rather to understand the periodic events of vegetation growth and development occurring at the pixel level.In the case of the savanna, it would entail studying the composite phenology of grasses and trees at the broad landscape.An extensive land surface phenological study undertaken by Guan et al. [1] on African savannas and woodlands revealed that African woodlands generally show early green-up before the onset of the rainy season and have a variable delayed senescence period after the rainy season.They also showed that the length of growing season had a positive linear relationship with tree cover fractions and that interannual variability in woodland phenological phases is smaller in comparison to grasses.The specific influence of varying tree cover fractions across the landscape on the spatio-temporal variations in land surface phenology in the semi-arid savanna of Southern Africa is not well understood.
Chidumayo [17] indicated that the ratio of tree:grass affects the spectral response (e.g., NDVI) of savanna pixels, and by implication should affect the phenological response of the pixels.For example, Fuller [22] showed that areas with relatively low tree canopy closure senesced earlier than tree dominated areas in Miombo woodlands.The main aim of this study was therefore to determine the influence of varying tree cover on the land surface phenology of the semi-arid savanna of Southern Africa, using the savanna landscape of South Africa as the study area.The influence of annual variation of rainfall on land surface phenology in the region was assessed.
Understanding the contribution of tree cover to land surface phenology could be relevant in the assessment of savanna ecosystem processes such as the spread of fire.The spread of fire in Southern African savanna is mainly determined by the amount of vegetation (herbaceous biomass and tree cover) and the state of the vegetation e.g., moisture content [23,24] and by implication, the phenological status of the vegetation.Bajocco et al. [25] concluded that the vegetation phenological status represents the main driver affecting fuel availability and moisture content and further argued that any investigation on fire behaviour over large areas requires the ability of mapping the spatio-temporal changes in vegetation phenology.Staver et al. [26] showed that fire spreads only at tree cover less than 40%.Therefore, the combined information on the contribution of trees to the land surface phenology of savannas and tree cover fractions could improve fire risk assessment in these systems.

Materials and Methods
The methodology involved: determining land surface phenological metric from leaf area index (LAI) time series data via a Gaussian model approach; -determining the influence of percentage tree cover on various land surface phenological metrics in two ways: (i) assessing whether significant differences exist between three tree cover classes for each land surface phenological metric; and (ii) assessing the strength of the relationship between each land surface phenological metric and percentage tree cover; and assessing the influence of rainfall on land surface phenology by Pearson's correlation between mean rainfall data with the land surface phenological metric.
This section starts by first describing the study area, and how the input data into the various analyses were derived i.e., the: (i) percentage tree cover for 100 randomly selected plots using high resolution Google Earth imagery; and (ii) LAI time series imagery by applying a regression model developed from synthetic data simulated using a canopy radiative transfer model on MODIS images.

Study Area
The study area consists of MODIS tile: h20v11, clipped to South African boundaries (Figure 1).Two main biomes dominate the MODIS tile, namely the savanna and grassland [9,27] (Figure 1a).The main focus though was on the savanna biome, which is in a hot semi-arid climatic zone, characterised by summer rainfall (September to April).The 2002-2015 precipitation data of five meteorological stations (Rustenburg, Marbel Hall, Moosriver, Lephalale, Bloemhof and Escelsior-benchmark) located in the region showed a mean annual rainfall of 820 mm.The rainfall data were obtained from the Agricultural Research Council (ARC) of South Africa.However, as can be observed in Figure 1b, mean annual rainfall varies across the region.(Schulze, 2007).One hundred random plots selected to analyse spatial and temporal patterns in land surface phenological metrics (c).

Tree Cover Percentage Estimation
We have assessed the quantity of woody plants in the system using tree canopy cover percentage [10].The tree cover estimates were determined from Google Earth high resolution imagery using similar methods described by Duhl et al. [28] and Ludwig et al. [29].One hundred fifty random polygons (2 by 2 pixels or 926 m by 926 m) were generated on the MODIS image scene, converted to kml files and overlaid on Google Earth high resolution image in the Google Earth software (Google Inc.) (Figure 1c).Only polygons that have not undergone changes in tree cover from 2001 to 2016 were retained for the study.This was achieved by examining historical images of each polygon.Polygons with no high resolution historical images were also not considered.Polygons that fell on plantation forest, crop fields and settlements were not considered and only images acquired during the summer months were used because woody cover in the savanna is known to vary with season [22].We stopped the exercise at the first 100 polygons.Three-band (red-green-blue) raster image scenes (JPEG format, Figure 2a) of the polygons were captured and imported into an image processing, ENVI software (ITT Visual information solutions).In ENVI, the scenes were classified into two classes (tree and non-tree) using a minimum distance supervised classification (see Figure 2 for an illustration of the process).The percentage tree cover per polygon was determined from the classified polygon as a ratio of the area of the tree-class to the total area of the polygon.(Schulze, 2007).One hundred random plots selected to analyse spatial and temporal patterns in land surface phenological metrics (c).

Tree Cover Percentage Estimation
We have assessed the quantity of woody plants in the system using tree canopy cover percentage [10].The tree cover estimates were determined from Google Earth high resolution imagery using similar methods described by Duhl et al. [28] and Ludwig et al. [29].One hundred fifty random polygons (2 by 2 pixels or 926 m by 926 m) were generated on the MODIS image scene, converted to kml files and overlaid on Google Earth high resolution image in the Google Earth software (Google Inc.) (Figure 1c).Only polygons that have not undergone changes in tree cover from 2001 to 2016 were retained for the study.This was achieved by examining historical images of each polygon.Polygons with no high resolution historical images were also not considered.Polygons that fell on plantation forest, crop fields and settlements were not considered and only images acquired during the summer months were used because woody cover in the savanna is known to vary with season [22].We stopped the exercise at the first 100 polygons.Three-band (red-green-blue) raster image scenes (JPEG format, Figure 2a) of the polygons were captured and imported into an image processing, ENVI software (ITT Visual information solutions).In ENVI, the scenes were classified into two classes (tree and non-tree) using a minimum distance supervised classification (see Figure 2 for an illustration of the process).The percentage tree cover per polygon was determined from the classified polygon as a ratio of the area of the tree-class to the total area of the polygon.

MODIS Imagery and Leaf Area Index (LAI) Time Series
Most studies have computed land surface phenological metrics from the normalised difference vegetation index (NDVI) or enhanced vegetation index (EVI) time series [30,31], but rarely from leaf area index (LAI) time series.The NDVI time series for African savannas are noted for their irregular temporal patterns [32].We also found the NDVI time series for the study area to be noisier than the LAI time series (e.g., Figure 3), even after applying a Savitzky-Golay smoothing filter (5 window filter, 5th order polynomial) to the data.The LAI time series for the MODIS tile (h20v11, MOD09A1 eight-day composite images at 463 m resolution) were established by using a regression model that was derived from synthetic data simulated from the PROSAIL radiative transfer model [33].PROSAIL is a combination of PROSPECT (a leaf optical properties model) and SAIL (a four-stream canopy radiative transfer model).First, synthetic canopy spectral reflectance data were simulated from the commonly used PROSAIL radiative transfer model in the 400-1060 nm optical range at 10 nm bandwidth, using MATLAB software.In the forward modelling mode, PROSAIL was parameterised using commonly used parameter ranges.Leaf (PROSPECT) parameters included: Chlorophyll content (0 to 90 µg/cm 2 ), equivalent water thickness (0.0005-0.055 cm), leaf dry matter content (1.25 × equivalent water thickness) and structural parameter (mean = 2.2 and standard deviation = 0.03, no dimension).Canopy (SAIL) parameters included: LAI (0 to 7.5 m 2 m −2 ), average leaf angle (25 to 80 degree), hotspot size (mean = 0.1 mm −1 , standard deviation = 0.01).The observation zenith angle, the solar zenith and azimuthal angle difference between the solar and view angles were fixed at 0, 0.78 rad and 0, respectively.The fraction of diffuse illumination was fixed at 0.1 independent of wavelengths.Two soil spectra (granite and dark gabbro soil) collected in the study area using the portable Analytical Spectral Device (ASD) were averaged and used to represent the soil spectral profile.Within the above parameter ranges, 1000 parameter sets were randomly generated and used to simulate the synthetic canopy reflectance data.
Two empirical models relating: (i) NDVI and LAI; and (ii) red and near-infrared bands, and LAI in a factorial regression analysis were explored using the synthetic data.The factorial regression model (Equation ( 1)) involving the red and NIR band was used to predict LAI on the MOD09A1 8day composite data for the period; January 2001 to December 2015, as this model yielded higher prediction accuracy (RMSE = 0.91) on a test dataset [33] compared to the regression model involving NDVI (RMSE = 1.33).Figure 4 shows a summer (Julian day 73) LAI map of South Africa.The LAI model adopted above showed a higher performance when compared to the MOD15 LAI product for the region [33].
where x and y are the red and NIR bands, respectively.The land surface phenological maps that were generated from the NDVI time series in our exploratory study contained far more speckles than those generated from the LAI time series.The LAI time series was therefore adopted in this study.Moreover, Archibald and Scholes [6] indicated

MODIS Imagery and Leaf Area Index (LAI) Time Series
Most studies have computed land surface phenological metrics from the normalised difference vegetation index (NDVI) or enhanced vegetation index (EVI) time series [30,31], but rarely from leaf area index (LAI) time series.The NDVI time series for African savannas are noted for their irregular temporal patterns [32].We also found the NDVI time series for the study area to be noisier than the LAI time series (e.g., Figure 3), even after applying a Savitzky-Golay smoothing filter (5 window filter, 5th order polynomial) to the data.The LAI time series for the MODIS tile (h20v11, MOD09A1 eight-day composite images at 463 m resolution) were established by using a regression model that was derived from synthetic data simulated from the PROSAIL radiative transfer model [33].PROSAIL is a combination of PROSPECT (a leaf optical properties model) and SAIL (a four-stream canopy radiative transfer model).First, synthetic canopy spectral reflectance data were simulated from the commonly used PROSAIL radiative transfer model in the 400-1060 nm optical range at 10 nm bandwidth, using MATLAB software.In the forward modelling mode, PROSAIL was parameterised using commonly used parameter ranges.Leaf (PROSPECT) parameters included: Chlorophyll content (0 to 90 µg/cm 2 ), equivalent water thickness (0.0005-0.055 cm), leaf dry matter content (1.25 × equivalent water thickness) and structural parameter (mean = 2.2 and standard deviation = 0.03, no dimension).Canopy (SAIL) parameters included: LAI (0 to 7.5 m 2 m −2 ), average leaf angle (25 to 80 degree), hotspot size (mean = 0.1 mm −1 , standard deviation = 0.01).The observation zenith angle, the solar zenith and azimuthal angle difference between the solar and view angles were fixed at 0, 0.78 rad and 0, respectively.The fraction of diffuse illumination was fixed at 0.1 independent of wavelengths.Two soil spectra (granite and dark gabbro soil) collected in the study area using the portable Analytical Spectral Device (ASD) were averaged and used to represent the soil spectral profile.Within the above parameter ranges, 1000 parameter sets were randomly generated and used to simulate the synthetic canopy reflectance data.
Two empirical models relating: (i) NDVI and LAI; and (ii) red and near-infrared bands, and LAI in a factorial regression analysis were explored using the synthetic data.The factorial regression model (Equation ( 1)) involving the red and NIR band was used to predict LAI on the MOD09A1 8-day composite data for the period; January 2001 to December 2015, as this model yielded higher prediction accuracy (RMSE = 0.91) on a test dataset [33] compared to the regression model involving NDVI (RMSE = 1.33).Figure 4 shows a summer (Julian day 73) LAI map of South Africa.The LAI model adopted above showed a higher performance when compared to the MOD15 LAI product for the region [33].LAI = 5.74 − 44.31x + 47.02x 2 − 3.43y + 4.72y 2 + 45.80xy where x and y are the red and NIR bands, respectively.The land surface phenological maps that were generated from the NDVI time series in our exploratory study contained far more speckles than those generated from the LAI time series.The LAI time series was therefore adopted in this study.Moreover, Archibald and Scholes [6] indicated that LAI would be preferable to NDVI in phenological studies, as it is one of the most accurate measures of the photosynthetic capacity of vegetation.
Remote Sens. 2017, 9, 689 6 of 19 that LAI would be preferable to NDVI in phenological studies, as it is one of the most accurate measures of the photosynthetic capacity of vegetation.

Determining Land Surface Phenological Metrics from LAI Time Series
Several methods have been used to extract phenological metrics from remotely sensed time series data including linear curve smoothing window and curve fitting methods [2,34,35].TIMESAT developed by Jonsson and Eklundh [34] is widely used software for deriving land surface phenological metrics.It is an open source software package that provides three smoothing functions that LAI would be preferable to NDVI in phenological studies, as it is one of the most accurate measures of the photosynthetic capacity of vegetation.

Determining Land Surface Phenological Metrics from LAI Time Series
Several methods have been used to extract phenological metrics from remotely sensed time series data including linear curve smoothing window and curve fitting methods [2,34,35].TIMESAT developed by Jonsson and Eklundh [34] is widely used software for deriving land surface phenological metrics.It is an open source software package that provides three smoothing functions

Determining Land Surface Phenological Metrics from LAI Time Series
Several methods have been used to extract phenological metrics from remotely sensed time series data including linear curve smoothing window and curve fitting methods [2,34,35].TIMESAT developed by Jonsson and Eklundh [34] is widely used software for deriving land surface phenological metrics.It is an open source software package that provides three smoothing functions to fit time series data for derivation of phenological metric including adaptive Savitzky-Golay filtering, asymmetric Gaussian and double logistic functions.Tan et al. [2] recommended the asymmetric Gaussian function against the other function because it is less sensitive to noise and missing data.In this study, we have used an inverted Gaussian function (f (t)) (Equation ( 2)), i.e., LAI as a function of date (t) of the growing season, that allows for the extraction of the most widely used phenological metrics as well as includes metrics such as day of maximum growing or senescing rates i.e., inflexion points on the growing (or green-up period) and senescing (senescence period) phases of the annual LAI curve (Figure 5).These include day of austral year at: (2) The above Gaussian function was therefore fitted to the entire LAI time series of each growing season which generally starts from 1 July (Julian Day 182, but 4 July in this study because of the MODIS 8-day interval) of one year to 30 June (Day 181; 26 of June in this study) of the following year.a = maximum LAI of the curve, b = left minimum of the LAI curve, c = day corresponding to the minimum LAI of the curve, σ = Gaussian variance and MGR = c + σ.There are no commonly agreed definitions of phenological metrics.For example, Wessel et al. [9] determined day of the year at SGS for the region used in this study as the date corresponding to 20% of season amplitude as measured from the left minimum of the curve.Many other definitions of day of the year at SGS and EGS are summarised in Groten and Ocatre [36].In this study, we have determined the day of the year at SGS for each growing season as the midpoint between the date of minimum LAI on the left side of the LAI time series and day of the year at MGR i.e., the period generally corresponding to the end of the lag phase of a growth curve.The day of the year at EGS and MSR were determined by using the same fitting procedure except applied in reverse order i.e., from 26 June backwards to 4 July of the previous year.Day of the year at PGS = (MGR + MSR)/2.A script was written in IDL (Interactive Data Language) to determine and map the phenological metrics for the study area.It involved making initial guesses of the model parameters; a, b, c and σ.The initial maps of the phenological metrics contained speckles which were smoothened out to a large degree by using a 3 by 3 window median filter.to fit time series data for derivation of phenological metric including adaptive Savitzky-Golay filtering, asymmetric Gaussian and double logistic functions.Tan et al. [2] recommended the asymmetric Gaussian function against the other function because it is less sensitive to noise and missing data.In this study, we have used an inverted Gaussian function (f(t)) (Equation ( 2)), i.e., LAI as a function of date (t) of the growing season, that allows for the extraction of the most widely used phenological metrics as well as includes metrics such as day of maximum growing or senescing rates i.e., inflexion points on the growing (or green-up period) and senescing (senescence period) phases of the annual LAI curve (Figure 5).These include day of austral year at: There are no commonly agreed definitions of phenological metrics.For example, Wessel et al. [9] determined day of the year at SGS for the region used in this study as the date corresponding to 20% of season amplitude as measured from the left minimum of the curve.Many other definitions of day of the year at SGS and EGS are summarised in Groten and Ocatre [36].In this study, we have determined the day of the year at SGS for each growing season as the midpoint between the date of minimum LAI on the left side of the LAI time series and day of the year at MGR i.e., the period generally corresponding to the end of the lag phase of a growth curve.The day of the year at EGS and MSR were determined by using the same fitting procedure except applied in reverse order i.e., from 26 June backwards to 4 July of the previous year.Day of the year at PGS = (MGR + MSR)/2.A script was written in IDL (Interactive Data Language) to determine and map the phenological metrics for the study area.It involved making initial guesses of the model parameters; a, b, c and σ.The initial maps of the phenological metrics contained speckles which were smoothened out to a large degree by using a 3 by 3 window median filter.

Data Analysis
The 100 tree cover polygons were grouped into three tree cover classes; Low (<20%, n = 44), Medium (20-40%, n = 22) and high (>40%, n = 34) to determine the influence of varying tree cover on the land surface phenological metrics of the semi-arid savanna of South Africa.We have determined the threshold for the tree cover classes based on expert opinion of the region (classification according to [37]), whereby areas with greater that 60% tree canopy cover (i.e., the canopies touch) in the savanna are referred to as forest, 40-60% as woodlands and less than 40% as open savanna.That is why we decided on three classes in the study i.e., evenly dividing 60% by three i.e., high (>40%), medium (20-40%) and low (0-20%).First, a median map (Figure 6) was generated for each metric from the 14 years considered in the study (2001/2002 to 2014/2015).Subsequently, the metric values of the various polygons were extracted from the median maps and the average for each polygon computed.Student's t-test was used to investigate if significant differences exist between the means of all pair-wise comparisons of tree cover classes i.e., low vs. medium, low vs. high and medium vs. high for all the phenological metrics.In addition, the interannual variabilities in the various metrics were assessed using the coefficient of variation (mean/Standard deviation × 100) for the 14 growing seasons.

Data Analysis
The 100 tree cover polygons were grouped into three tree cover classes; Low (<20%, n = 44), Medium (20-40%, n = 22) and high (>40%, n = 34) to determine the influence of varying tree cover on the land surface phenological metrics of the semi-arid savanna of South Africa.We have determined the threshold for the tree cover classes based on expert opinion of the region (classification according to [37]), whereby areas with greater that 60% tree canopy cover (i.e., the canopies touch) in the savanna are referred to as forest, 40-60% as woodlands and less than 40% as open savanna.That is why we decided on three classes in the study i.e., evenly dividing 60% by three i.e., high (>40%), medium (20-40%) and low (0-20%).First, a median map (Figure 6) was generated for each metric from the 14 years considered in the study (2001/2002 to 2014/2015).Subsequently, the metric values of the various polygons were extracted from the median maps and the average for each polygon computed.Student's t-test was used to investigate if significant differences exist between the means of all pair-wise comparisons of tree cover classes i.e., low vs. medium, low vs. high and medium vs. high for all the phenological metrics.In addition, the interannual variabilities in the various metrics were assessed using the coefficient of variation (mean/Standard deviation × 100) for the 14 growing seasons.Precipitation data of the five meteorological stations (Rustenburg, Marbel Hall, Moosriver, Lephalale, Bloemhof and Escelsior-benchmark) were obtained from the Agricultural Research Council (ARC) of South Africa for the period 2002 to 2015 (13 years) and used to investigate the impact of rainfall on the land surface phenological metrics.The start and end of rainy season have been used in previous studies to determine the impact of rainfall on land surface phenology e.g., Suepa et al. (2016) [3].However, there is no consensus on the definitions of the start and end of the rainy season [38].Furthermore, increasing rainfall variability at the start of the rainy season in Southern Africa provides an extra challenge in defining the start of the rainy season.We have therefore used the cumulative monthly mean rainfall in this study as the explanatory variable for variations in land surface phenology.In this respect, the rainfall data for each growing season or Austral year were processed into several categories (Table 1), consisting of cumulative monthly mean rainfall for the green-up phase of the vegetation i.e., July to December and senescence phase in a reverse order i.e., June-January.The relationships between the various rainfall categories and the land surface phenological metrics were subsequently assessed using the Pearson's correlation coefficient (r).

Spatial Patterns of Land Surface Phenological Metrics and Differences between Tree Cover (%) Classes
The median maps of the land surface phenological metrics show varying spatial patterns between the metrics on the growing (day of the year at SGS, MGR or PGS) and senescing (day of the year at MSR and EGS) phases of the seasonal LAI curve (Figure 6).The spatial patterns of the metrics on the growing phase closely follow the precipitation pattern of the study area, with early day of the year at SGS, MGR or PGS regions associated with the regions of high annual mean precipitation (compare Figures 1 and 6).While the maps of the metrics on the senescing phase show similar patterns to the summer LAI variability across the study area (compare Figures 3 and 6).The following general timings of the growing season can be associated with the study area using analysis of the 100 randomly selected polygons and following our definition of the phenological metrics (Table 2) The vegetated landscape starts to green-up (day of the year at SGS) at about August 31 ± 2 days (95% CI). - The vegetated landscape reaches peak production (PGS) at about January 20 ± 3 days (95% CI). - The growing season ends (EGS) at about May 11 ± 1 day (95% CI).

Differences in Land Surface Phenological Metrics between Tree Cover Classes
The 100 randomly selected polygons grouped into three percentage tree cover classes were used in this analysis.First, the differences in the means between tree cover classes were not statistically significant (t-test p > 0.05) for each of the three growing phase phenological metrics (day of the year at SGS, MGR and PGS) (Table 3, Figure 7), though with higher means for the low tree cover class for day of the year at SGS and MGR in comparison to the medium and high tree cover classes.Secondly, the variability in the low tree cover class for the growing phase metrics was considerably higher (average CV = 19.4%,p variances < 0.05) than in the medium (CV = 11.7%) or high (CV = 9.7%) tree cover class.In fact, the range of values of the low tree class includes the entire range of values for the medium and high tree covers classes.On the other hand, the three tree cover classes were clearly distinguishable (t-test p < 0.01) using the senescing phase metrics (day of the year at MSR and EGS).The coefficients of variation of the senescing phase metrics (average CV = 2%) were far lower than those of the growing phase metrics (average CV = 13.6%).
To further understand the influence of tree cover on the land surface phenological metrics, we conducted linear regression analyses between tree cover and the various land surface phenological metrics (Figure 8).Tree cover (%) showed a slightly non-significant negative correlation with day of the year at SGS, MGR, or PGS (p > 0.05), but a highly positive correlation with day of the year at MSR or EGS (p < 0.01).In fact, the relationship between tree cover and day of the year at MSR or EGS is non-linear, with a sharp rise in day of the year at MSR or EGS between 0% and about 10% tree cover.A Gaussian model was fitted (R 2 = 0.40) to the scatter plot between day of the year at MSR and percentage tree cover (Figure 9) and the model parameters were used to produce a tree cover map for the study area (Figure 10).Table 3. Differences in the means and variances between tree cover (%) classes using t and f tests.Table 3. Differences in the means and variances between tree cover (%) classes using t and f tests.by almost 14 days, while the landscape showed an earlier day of the year at SGS in 2007 with an above average July-September mean rainfall.

Discussion
The main aim of this study was to determine the influence of varying tree cover on land surface phenology of the semi-arid savanna in South Africa.In other words, what is the impact of variation in tree cover percentage on the phenological response of mixed remotely sensed pixel of the semiarid savanna?
An important finding of this study is the differential response of the land surface phenology to tree cover during the growing and senescing phases of the vegetated landscape.The influence of tree cover on land surface phenology is far stronger on the senescing phase metrics than on the growing phase metrics, e.g., explaining about 40% of the variance (p < 0.05) in the day of the year at MSR, i.e., point of inflexion on the senescing phase of the austral year.This study thus suggests an alternative method of assessing tree cover fraction in the semi-arid savanna i.e., using senescence period land surface phenological metrics.The varying influence of tree cover on the land surface phenology is most probably dependent on the vegetation response to climatic variables such as rainfall and temperature.Rainfall was shown to have a strong influence on land surface phenological metrics

Discussion
The main aim of this study was to determine the influence of varying tree cover on land surface phenology of the semi-arid savanna in South Africa.In other words, what is the impact of variation in tree cover percentage on the phenological response of mixed remotely sensed pixel of the semi-arid savanna?
An important finding of this study is the differential response of the land surface phenology to tree cover during the growing and senescing phases of the vegetated landscape.The influence of tree cover on land surface phenology is far stronger on the senescing phase metrics than on the growing phase metrics, e.g., explaining about 40% of the variance (p < 0.05) in the day of the year at MSR, i.e., point of inflexion on the senescing phase of the austral year.This study thus suggests an alternative method of assessing tree cover fraction in the semi-arid savanna i.e., using senescence period land surface phenological metrics.The varying influence of tree cover on the land surface phenology is most probably dependent on the vegetation response to climatic variables such as rainfall and temperature.Rainfall was shown to have a strong influence on land surface phenological metrics derived from the growing phase of the vegetation, thus overshadowing the impact of tree cover on the phenological response.Moreover, tree green-up has been shown to vary amongst savanna tree species with some tree species greening up before and others after the onset of the rainy season [6].Thus, the variability in rainfall across the region and the differential behaviour of tree species in the region in relation to the onset of the rainy season might explain the low relationship between tree cover and start of the growing season.Temperature as opposed to rainfall might be the dominant factor controlling senescence of the vegetated landscape.
There were higher variabilities across the region (i.e., using the 100 randomly selected plots), for growing phase metrics when compared to the senescing phase metrics.Again, we believe that the spatial variability of rainfall in the region accounts for the high variabilities in the growing phase metrics.In fact, the low tree cover class showed the highest variabilities in the growing phase metrics.The low tree areas are most probably dominated by grasses and it has been widely established that grass phenology is largely dependent on rainfall [1,17].Therefore, as asserted by Chidumayo [17], spatial responses of the start of the growing season can be explained by differential responses of phenology to rainfall among grasses and trees.Land surface phenological metrics derived from the growing phase of the annual LAI time series showed higher interannual variability in comparison to the senescence period metrics.Interannual variability of rainfall also explains the interannual variability of the growing phase metrics which were shown to be more sensitive to variation in rainfall.Contrary to expectations, we did not find any differences in the interannual variability between the tree cover classes at any phenological stage.Tree phenology was found to be more stable between years when compared to the phenology of savanna grasses in the Southern part of the Kruger National Park, which is within the study area [6].Guan et al. [1] also determined that interannual variability in woodland phenological phases is smaller in comparison to African grasslands.The impact of rainfall on the interannual variability of land surface phenological metrics in the region may be site dependent.This however remains to be established.It is also noteworthy to mention that several challenges in retrieving the land surface phenological metrics could have affected the results.These challenges include cloud covered areas which cause distortion in the seasonal profile, no clear season LAI profiles in certain areas e.g., in a close canopy environment or very dry areas or several amplitudes in some areas which may results from subtle green-up events.Furthermore, the LAI time series maps were generated by applying a regression model derived from the red and NIR reflectance bands that were simulated using the PROSAIL radiative transfer model (a 1-D model).We however recommend further investigation into how the results might be different if LAI is generated using other radiative transfer models and inversion techniques.
There was little inter-annual variability in the senescing phase metrics.We speculate that temperature plays a more important role during this later phase than rainfall.There may be less spatial and temporal variation in temperature during the senescing phase thereby impacting less on the senescing dates of the various tree cover classes, reason for which the three tree cover classes were clearly distinguishable during senescence.Chidumayo [17] concluded that the interaction between minimum and maximum temperature was the most important determinant of savanna phenology in Southern Africa.Chidumayo [39] further determined that the interactions between temperature and rainfall explained 60-98% of the variations in diameter at breast height increment in indigenous trees in Zambia.Guan et al. [1] however concluded that potential factors controlling the phenology of African woodlands are not yet clearly understood and suggested further exploration of this subject.The complexity in understanding land surface phenology of savanna woodlands is further compounded by the mixed nature of remotely sensed pixels.
We argue that understanding the contribution of tree cover to land surface phenology could be crucial to understanding ecological processes such as the spread of fire in savanna ecosystems.Fire occurs in the African savanna system mostly at the end of the growing season and its spread is mainly determined by the dry herbaceous biomass and vegetation state e.g., moisture content of the vegetation or the phenological status.We now know that that spatial variability of day of the year at end of growing season is mainly determined by the fraction of tree cover in the savanna system of South Africa.Thus, the map of the spatial variability of day of the year at end of growing season could be important in assessing the spatial variability in fire spread with time.Furthermore, the low interannual variability in the day of the year at end of growing season implies that beyond this day for grass dominated areas (at about Day 300 of the austral year), any further spatial variation in moisture content of the vegetation with time could most probably be attributed to variation in tree cover fraction.We recommend the above assertion for further investigation.

Conclusions
The main results of the study can be summarised as follows:

•
Areas with tree cover greater than 20% showed an earlier day of the year at SGS or MGR in comparison to areas with tree cover of less than 20%, although the difference was statistically not significant at p < 0.05 due to the wide variation in day of the year at SGS or MGR among the low tree cover polygons located across the study region.

•
There was a significant positive trend in the phenological metrics computed from the senescing phase, e.g., day of the year at MSR or EGS increases with increasing tree cover.

•
Phenological metrics computed from the growing phase of the annual LAI time series showed higher inter-annual variability in comparison to the senescing phase metrics.

•
Rainfall data for the region explained inter-annual variations in land surface phenological metrics during the green-up period of the vegetation, but not for the senescence period.
In conclusion, land surface phenological metrics derived from the green-up and senescence periods of the semi-arid savanna of Southern Africa exhibit differential response to tree cover variations.Variation in land surface phenology in the green-up period is partly dependent on rainfall and mean annual precipitation varies across the region [40], while tree cover strongly influences senescence period phenological response.The combined information on the contribution of trees to the land surface phenology of savanna and tree cover fractions could be relevant in models of fire risk assessment across the region with time since the spread of fire is determined by fuel load and the state of the vegetation e.g., phenology.

Figure 1 .
Figure 1.Biomes of South Africa (a) (Mucina and Rutherford, 2006); and mean annual precipitation (b) (Schulze, 2007).One hundred random plots selected to analyse spatial and temporal patterns in land surface phenological metrics (c).

Figure 1 .
Figure 1.Biomes of South Africa (a) (Mucina and Rutherford, 2006); and mean annual precipitation (b) (Schulze, 2007).One hundred random plots selected to analyse spatial and temporal patterns in land surface phenological metrics (c).

Figure 2 .
Figure 2. Tree cover (%) estimation using classification of high resolution Google Earth imagery (a) into tree and background classes (b).The polygon boundary is indicated in black on image (a).

Figure 2 .
Figure 2. Tree cover (%) estimation using classification of high resolution Google Earth imagery (a) into tree and background classes (b).The polygon boundary is indicated in black on image (a).

Figure 3 .
Figure 3. Sample time series of two spectral indices extracted from MODIS imagery of the study area, (a) the Normalised Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI), and (b) leaf area index (LAI) derived from the red and near-infrared bands.

Figure 4 .
Figure 4. Leaf area index map of South Africa (Julian day 73) generated by inversion of PROSAIL radiative transfer model on MODIS imagery.

Figure 3 .
Figure 3. Sample time series of two spectral indices extracted from MODIS imagery of the study area, (a) the Normalised Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI), and (b) leaf area index (LAI) derived from the red and near-infrared bands.

Figure 3 .
Figure 3. Sample time series of two spectral indices extracted from MODIS imagery of the study area, (a) the Normalised Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI), and (b) leaf area index (LAI) derived from the red and near-infrared bands.

Figure 4 .
Figure 4. Leaf area index map of South Africa (Julian day 73) generated by inversion of PROSAIL radiative transfer model on MODIS imagery.

Figure 4 .
Figure 4. Leaf area index map of South Africa (Julian day 73) generated by inversion of PROSAIL radiative transfer model on MODIS imagery.
start of growing season (SGS) • end of growing season (EGS) • peak or maximum growth (PGS) • maximum growing rate (MGR) i.e., the inflexion point on the growing phase • maximum senescing rate (MSR) i.e., inflexion point on the senescing phase (2) The above Gaussian function was therefore fitted to the entire LAI time series of each growing season which generally starts from 1 July (Julian Day 182, but 4 July in this study because of the MODIS 8-day interval) of one year to 30 June (Day 181; 26 of June in this study) of the following year.a = maximum LAI of the curve, b = left minimum of the LAI curve, c = day corresponding to the minimum LAI of the curve, σ = Gaussian variance and MGR = c + σ.

Figure 5 .
Figure 5. Land surface phenological metrics resulting from fitting an inverted Gaussian model to the LAI austral year curve on the green-up and senescence periods of the vegetated landscape.SGS, MGR, PGS, MSR and EGS denote dates corresponding to start, maximum growing rate, peak growth, maximum senescing rate and end of growing season, respectively.

Figure 5 .
Figure 5. Land surface phenological metrics resulting from fitting an inverted Gaussian model to the LAI austral year curve on the green-up and senescence periods of the vegetated landscape.SGS, MGR, PGS, MSR and EGS denote dates corresponding to start, maximum growing rate, peak growth, maximum senescing rate and end of growing season, respectively.

Figure 6 .
Figure 6.Fourteen-year median maps of the various phenological metrics: (a-e) the day of the austral year at start (SGS), maximum growing rate (MGR), peak growth (PGS), maximum senescing rate (MSR) and end of growing season (EGS), respectively.The austral year used in the study ranges from 4 July (Julian Day 185) of one year to 26 June (Day 177) of the following year.

Figure 6 .
Figure 6.Fourteen-year median maps of the various phenological metrics: (a-e) the day of the austral year at start (SGS), maximum growing rate (MGR), peak growth (PGS), maximum senescing rate (MSR) and end of growing season (EGS), respectively.The austral year used in the study ranges from 4 July (Julian Day 185) of one year to 26 June (Day 177) of the following year.

Figure 7 .
Figure 7. Differences in land surface phenology between three tree cover (%) classes (low, medium and high) for day of the year at: (a) start of growing season; (b) maximum growing rate; (c) peak growth; (d) maximum senescing rate; and (e) end of growing season.The Austral used in the study ranges from 4 July (Julian Day 185) of one year to 26 June (Day 177) of the following year.

Figure 7 .
Figure 7. Differences in land surface phenology between three tree cover (%) classes (low, medium and high) for day of the year at: (a) start of growing season; (b) maximum growing rate; (c) peak growth; (d) maximum senescing rate; and (e) end of growing season.The Austral used in the study ranges from 4 July (Julian Day 185) of one year to 26 June (Day 177) of the following year.

Figure 8 .
Figure 8. Linear regression analyses between tree cover (%) and landscape phenological metrics for day of the year at: (a) start of growing season; (b) maximum growing rate; (c) peak production; (d) maximum senescing rate; and (e) end of growing season.The growing season used in the study ranges from 4 July (Julian Day 185) of one year to 26 June (Day 177) of the following year.

Figure 8 .
Figure 8. Linear regression analyses between tree cover (%) and landscape phenological metrics for day of the year at: (a) start of growing season; (b) maximum growing rate; (c) peak production; (d) maximum senescing rate; and (e) end of growing season.The growing season used in the study ranges from 4 July (Julian Day 185) of one year to 26 June (Day 177) of the following year.

Figure 9 .
Figure 9. Model used to map tree cover (%) in the study area.

Figure 10 .
Figure 10.Tree cover (%) derived from the day of the austral year at maximum senescing rate (a land surface phenological metric).

Figure 9 .
Figure 9. Model used to map tree cover (%) in the study area.

Figure 9 .
Figure 9. Model used to map tree cover (%) in the study area.

Figure 10 .
Figure 10.Tree cover (%) derived from the day of the austral year at maximum senescing rate (a land surface phenological metric).

Figure 10 .
Figure 10.Tree cover (%) derived from the day of the austral year at maximum senescing rate (a land surface phenological metric).

Figure 12 .
Figure 12.Yearly deviation of July-September mean rainfall (5 meteorological stations) from the 13 year (2002-2014) mean (a); and yearly deviation of start of growing season date for 100 randomly selected polygons in the study area (b).

Figure 12 .
Figure 12.Yearly deviation of July-September mean rainfall (5 meteorological stations) from the 13 year (2002-2014) mean (a); and yearly deviation of start of growing season date for 100 randomly selected polygons in the study area (b).

Table 1 .
Rainfall data of five meteorological stations in the study area (Rustenburg, Marbel Hall Moosriver, Lephalale, Bloemhof and Escelsior-benchmark), showing mean (mm) rainfall for the green-up and senescence periods of the vegetated landscape.

Phenological Metrics t-Value for Differences in Means p for Differences in Means F-Ratio for Differences in Variances p for Differences Variances
Bold p values are significant at p < 0.05.

Table 4 .
Pearson's correlations (r) between mean rainfall of five meteorological stations in the study area (Rustenburg, Marbel Hall Moosriver, Lephalale, Bloemhof and Escelsior-benchmark) and green-up or senescence period phenological metrics for 100 randomly selected plots in the study area.

Table 4 .
Pearson's correlations (r) between mean rainfall of five meteorological stations in the study area (Rustenburg, Marbel Hall Moosriver, Lephalale, Bloemhof and Escelsior-benchmark) and greenup or senescence period phenological metrics for 100 randomly selected plots in the study area.