Deriving seasonal dynamics in ecosystem properties of semi-arid savanna grasslands from in situ-based hyperspectral reflectance

. This paper investigates how hyperspectral reﬂectance (between 350 and 1800 nm) can be used to infer ecosystem properties for a semi-arid savanna grassland in West Africa using a unique in situ-based multi-angular data set of hemispherical conical reﬂectance factor (HCRF) measurements. Relationships between seasonal dynamics in hyperspectral HCRF and ecosystem properties (biomass, gross primary productivity (GPP), light use efﬁciency (LUE), and fraction of photosynthetically active radiation absorbed by vegetation (FAPAR)) were analysed. HCRF data ( ρ ) were used to study the relationship between normalised difference spectral indices (NDSIs) and the measured ecosystem properties. Finally, the effects of variable sun sensor viewing geometry on different NDSI wavelength combinations were analysed. The wavelengths with the strongest correlation to seasonal dynamics in ecosystem properties were shortwave infrared (biomass), the peak absorption band for chlorophyll a and b (at 682 nm) (GPP), the oxygen A band at 761 nm used for estimating chlorophyll ﬂuorescence (GPP and LUE), and blue wavelengths ( ρ 412 ) (FAPAR). The NDSI with the strongest correlation to (i) biomass combined red-edge HCRF ( ρ 705 ) with green HCRF ( ρ 587 ), (ii) GPP combined wavelengths at the peak of green reﬂection ( ρ 518 , ρ 556 ), (iii) LUE combined red ( ρ 688 ) with blue HCRF ( ρ 436 ), and (iv) FAPAR combined blue ( ρ 399 ) and near-infrared ( ρ 1295 ) wavelengths. NDSIs combining near infrared and shortwave infrared were strongly affected by solar zenith angles and sensor viewing geometry, as were many combinations of visible wavelengths. This study provides analyses based upon novel multi-angular hyperspectral data for validation of Earth-observation-based properties of semi-arid ecosystems, as well as insights for designing spectral characteristics of future sensors for ecosystem monitoring.


Introduction
Hyperspectral measurements of the Earth's surface provide relevant information for many ecological applications. An important tool for spatial extrapolation of ecosystem functions is to study how spectral properties are related to in situ measured ecosystem properties. These relationships found the basis for upscaling using Earth observation (EO) data. Continuous in situ measurements of hyperspectral reflectance in combination with ecosystem properties are thereby essential for improving our understanding of the functioning of the observed ecosystems. Strong relationships have, for example, been found between information in the reflectance spectrum and ecosystem properties such as leaf area index (LAI), fraction of photosynthetically active radiation (PAR) absorbed by the vegetation (FAPAR), light use efficiency (LUE), biomass, vegetation primary productivity, vegetation water content, and nitrogen and chlorophyll content (e.g. Thenkabail et al., 2012;Tagesson et al., 2009; Published by Copernicus Publications on behalf of the European Geosciences Union.

4622
T. Tagesson et al.: Deriving seasonal dynamics in ecosystem properties of semi-arid savanna grasslands Gower et al., 1999;Sjöström et al., 2009;Sims and Gamon, 2003). In situ observations of spectral reflectance are also important for parameterisation and validation of canopy reflectance models, as well as space-and airborne products (Coburn and Peddle, 2006).
Very few sites across the world exist with an instrumental setup designed for multi-angular continuous hyperspectral measurements. Leuning et al. (2006) present a system mounted in a 70 m tower above an evergreen eucalyptus forest in New South Wales, Australia, which measures spectral hemispherical conical reflectance factors (HCRFs) 1 hourly throughout the year between 300 and 1150 nm at four azimuth angles. Hilker et al. (2007Hilker et al. ( , 2010 describe an automated multi-angular spectro-radiometer for estimation of canopy HCRF (AMSPEC) mounted on a tower above a coniferous forest in Canada. Spectral HCRF is sampled between 350 and 1200 nm year-round under different viewing and sun angle conditions, achieved by collection of data in a near 360 • view around the tower with adjustable viewing zenith angles. Even though in situ measurements of multiangular hyperspectral HCRF are fundamental for the EO research community, such data sets are still rare and, at the present state, they do not cover different biomes at the global scale (Huber et al., 2014).
There are many methods for analysing relationships between hyperspectral reflectance and ecosystem properties, such as multivariate methods, derivative techniques, and radiative transfer modelling (Bowyer and Danson, 2004;Ceccato et al., 2002;Danson et al., 1992;Roberto et al., 2012). Still, due to its simplicity, the combination of reflectance into vegetation indices is the major method for upscaling using EO data. By far, the most commonly applied vegetation indices are those based on band ratios, e.g. the normalised difference vegetation index (NDVI), which is calculated by dividing the difference in the near-infrared (NIR) and red wavelength bands by the sum of the NIR and red bands (Tucker, 1979;Rouse et al., 1974). The NIR radiance is strongly scattered by the air-water interfaces between the cells, whereas red radiance is absorbed by chlorophyll and its accessory pigments (Gates et al., 1965). The normalisation with the sum in the denominator is a means to reduce the effects of solar zenith angle, sensor viewing geometry, and atmospheric errors as well as enhancing the signal of the observed target (e.g. Qi et al., 1994;Inoue et al., 2008).
Wavelength specific spectral reflectance is known to be related to leaf characteristics such as chlorophyll concentra-1 Different reflectance terminologies have been used to inform on spectral measurements in the field by the remote sensing community leading to suggestions to the proper use of the terminology (Martonchik et al., 2000). All field spectro-radiometers measure HCRF (hemispherical conical reflectance) if the field of view (FOV) of the sensor is larger than 3 • (Milton et al., 2009). HCRF is therefore used throughout this paper to support the correct inference and usage of reflectance products (Schaepman-Strub et al., 2006;Milton et al., 2009). tion, dry matter content, internal structure parameters, and equivalent water thickness (Ceccato et al., 2002). Hyperspectral reflectance data can be combined into a matrix of normalised difference spectral indices (NDSIs), following the NDVI rationing approach. Correlating the NDSI with ecosystem properties provides a way to gain an improved empirically based understanding of the relationship between information in the reflectance spectrum with ground surface properties (e.g. Inoue et al., 2008). Several studies have analysed relationships between hyperspectral HCRF, NDSI, and ecosystem properties (e.g. Thenkabail et al., 2000;Cho et al., 2007;Psomas et al., 2011;Inoue et al., 2008;Gamon et al., 1992;Feret et al., 2008;Thenkabail et al., 2012). Still, it is extremely important to examine these relationships for different ecosystems across the Earth and investigate their applicability for different environmental conditions and under different effects of biotic and abiotic stresses.
A strong correlation between an NDSI and an ecosystem property does not necessarily indicate that the NDSI is a good indicator of vegetation conditions to be applied to EO systems. Visible, NIR, and shortwave infrared (SWIR) have different sensitivity to variations in solar zenith angles, stand structure, health status of the vegetation, vegetation and soil water content, direct/diffuse radiation ratio, and sensor viewing geometry. The influence of sun-sensor geometry on the reflected signal has been studied using radiative transfer models and airborne (e.g. AirMISR) as well as satellitebased data from instruments such as CHRIS-PROBA, MISR, or POLDER (Huber et al., 2010;Maignan et al., 2004;Javier García-Haro et al., 2006;Jacquemoud et al., 2009;Verhoef and Bach, 2007;Laurent et al., 2011). However, effects of variable sun angles and sensor viewing geometries are not well documented in situ for different plant functional types of natural ecosystems except for some individual controlled experiments (Hilker et al., 2008;Sandmeier et al., 1998;Schopfer et al., 2008). Improved knowledge regarding the influence from sun-sensor variability on different NDSI combinations is thereby essential for validating the applicability of an NDSI for EO upscaling purposes.
The Dahra field site in Senegal, West Africa, was established in 2002 as an in situ research site to improve our knowledge regarding properties of semi-arid savanna ecosystems and their responses to climatic and environmental changes (Tagesson et al., 2015b). A strong focus of this instrumental setup is to gain insight into the relationships between ground surface reflectance and savanna ecosystem properties for EO upscaling purposes. This paper presents a unique in situ data set of seasonal dynamics in hyperspectral HCRF and demonstrates how it can be used to describe the seasonal dynamics in ecosystem properties of semi-arid savanna ecosystems. The objectives are threefold: (i) to quantify the relationship between seasonal dynamics of in situ hyperspectral HCRF between 350 and 1800 nm and ecosystem properties (biomass, gross primary productivity (GPP), LUE, and FAPAR), (ii) to quantify the relationship between NDSI Table 1. Information about the instrumental setup for the measured environmental variables. HCRF is hemispherical conical reflectance factor, GPP is gross primary productivity, LUE is light use efficiency, and FAPAR is fraction of photosynthetically active radiation absorbed by the vegetation. Min and Max are minimum and maximum values measured, respectively, DW is dry weight, C is carbon, and MJ is megajoules. The year started is the first year with measurements. Time is in UTC. For more information about the instrumental setup, see Tagesson et al. (2015b).

Variable
Year with different wavelength combinations (350 to 1800 nm) and the measured ecosystem properties, and (iii) to analyse and quantify effects of variable sun angles and sensor viewing geometries on different NDSI combinations.

Site description
All measurements used for the present study were conducted at the Dahra field site in the Sahelian ecoclimatic zone northeast of the town Dahra in the semi-arid central part of Senegal (15 • 24 10 N, 15 • 25 56 W) during 2011 and 2012 (Fig. 1). Rainfall is sparse in the region, with a mean annual sum of 416 mm . More than 95 % of the rain falls between July and October, with August being the wettest month. The mean annual air temperature is 29 • C (1951-2003); May is the warmest month and January the coldest, with mean monthly temperatures of 32 and 25 • C, respectively. The Dahra site has a short growing season (∼ 3 months) following the rainy season, with leaf area index generally ranging between 0 and 2 (Fensholt et al., 2004). South-western winds dominate during the rainy season and north-eastern winds during the dry season. The area is dominated by annual grasses (e.g. Schoenefeldia gracilis, Digitaria gayana, Dactyloctenium aegypticum, Aristida mutabilis, and Cenchrus biflorus) (Mbow et al., 2013), and trees and shrubs (e.g. Acacia senegal and Balanites aegyptiaca) are relatively sparse (∼ 3 % of the land cover) (Rasmussen et al., 2011). The average tree height was 5.2 m and the peak height of the herbaceous layer was 0.7 m (Tagesson et al., 2015b). A thorough description of the Dahra field site is given in Tagesson et al. (2015b).

Meteorological and vegetation variables
A range of meteorological variables have been measured from a tower at the Dahra field site for more than 10 years: air temperature ( • C) and relative humidity (%) were mea-sured at 2 m height; soil temperature ( • C) and soil moisture (volumetric water content (m 3 m −3 × 100) (%)) were collected at 0.05 m depths; rainfall (mm) was measured at 2 m height; incoming ( inc ) and reflected ( ref ) PAR (µmol m −2 s −1 ) was measured at 10.5 m height; and PAR transmitted through the vegetation (PAR transmit ) was measured at six plots at ∼ 0.01 m height (Table 1) (Tagesson et al., 2015b). The PAR transmit was measured within a distance of 7 m from the tower. PAR absorbed by the vegetation (APAR) was estimated by where α soil is the PAR albedo of the soil, which was measured as 0.20 (Tagesson et al., 2015b). FAPAR was estimated by dividing APAR by PAR inc (Tagesson et al., 2015b). All sensors were connected to a CR-1000 logger in combination with a multiplexer (Campbell Scientific Inc., North Logan, USA). Data were sampled every 30 s and stored as 15 min averages (sum for rainfall).
The total above-ground green biomass (g m −2 ) of the herbaceous vegetation was sampled approximately every 10 days during the growing seasons 2011 and 2012 at twentyeight 1 m 2 plots located along two ∼ 1060 m long diagonal transects ( Fig. 1f) (Mbow et al., 2013). The method applied was destructive, so even though the same transects were used for each sampling date, the plots were never positioned at exactly the same location. The study area is flat and characterised by homogenous grassland savanna and the conditions in these sample plots are generally found to be representative for the conditions in the entire measurement area (Fensholt et al., 2006). All above-ground green herbaceous vegetation matter was collected and weighed in the field to get the fresh weight. The dry matter (DW) was estimated by oven-drying the green biomass. For a thorough description regarding the biomass sampling we refer the reader to Mbow et al. (2013). showing the location of the meteorological tower, the EC tower, the biomass sampling plots, and the footprint of the EC measurements. The EC footprint area is the median 70 % cumulative flux distance from the eddy covariance tower. The photos of the EC tower and its footprint are taken during the rainy season, whereas the photo of the meteorological tower shows the late dry season.

Estimates of gross primary productivity and light use efficiency
Net ecosystem exchange of CO 2 (NEE) (µmol CO 2 m −2 s −1 ) was measured with an eddy covariance system, consisting of an open-path infrared gas analyser (LI-7500, LI-COR Inc., Lincoln, USA) and a three-axis sonic anemometer (Gill instruments, Hampshire, UK) from 18 July 2011 until 31 December 2012 (Table 1). The sensors were mounted 9 m above the ground on a tower (placed 50 m south of the tower including the meteorological and spectroradiometric sensors) ( Fig. 1b and f). Data were sampled at 20 Hz rate. The postprocessing was done with the EddyPro 4.2.1 software (LI-COR Biosciences, 2012), and statistics were calculated for 30 min periods. The post-processing includes 2-D coordinate rotation (Wilczak et al., 2001), time lag removal between anemometer and gas analyser by covariance maximisation (Fan et al., 1990), despiking (Vickers and Mahrt, 1997) (plausibility range: window average ±3.5 SD), linear detrending (Moncrieff et al., 2004), and compensation for density fluctuations (Webb et al., 1980). Fluxes were also corrected for high-pass (Moncrieff et al., 1997) and low-pass filtering ef- fects (Moncrieff et al., 2004). The data were filtered for steady state and fully developed turbulent conditions, following Foken et al. (2004), and according to statistical tests as recommended by Vickers and Mahrt (1997). Flux measurements from periods of heavy rainfall were also removed. For a thorough description of the post processing of the raw eddy covariance data, see Tagesson et al. (2015a). A possible source of error in a comparison between ECbased variables and spectral HCRF is the difference in footprint/instantaneous field of view (IFOV) between the sensors. The IFOV of the spectroradiometer setup contains only soil and herbaceous vegetation. The footprint of the EC tower was estimated using a model based on measurement height, surface roughness, and atmospheric stability (Hsieh et al., 2000). The median point of maximum contribution is at 69 m, and the median 70 % cumulative flux distance is at 388 m from the tower. The footprint of the EC tower contains semi-arid savanna grassland with ∼ 3 % tree coverage, and the EC data are thereby affected by both woody and herbaceous vegetation ( Fig. 1a and f). But given the low tree coverage, and the dominant influence of herbaceous vegetation on the seasonal dynamics in CO 2 fluxes, we still consider it reasonable to compare EC fluxes with seasonal dynamics in spectral HCRF of the herbaceous vegetation.

Biogeosciences
The daytime NEE was partitioned into GPP and ecosystem respiration using the Mitscherlich light-response function against PAR inc (Falge et al., 2001). A 7-day moving window with 1-day time steps was used when fitting the functions. By subtracting dark respiration (R d ) from the lightresponse function, it was forced through 0, and GPP was estimated: where F csat is the CO 2 uptake at light saturation (µmol CO 2 m −2 s −1 ) and α is the quantum efficiency or the initial slope of the light-response curve (µmol CO 2 (µmol photons) −1 ) (Falge et al., 2001). Vapour pressure deficit (VPD) limits GPP and to account for this effect, the F csat parameter was set as an exponentially decreasing function: where VPD 0 is 10 hPa following the method by Lasslop et al. (2010). Gaps in GPP less than or equal to 3 days were filled using three different methods: (i) gaps shorter than 2 h were filled using linear interpolation, (ii) daytime gaps were filled by using the light-response function for the 7-day moving windows, and (iii) remaining gaps were filled by using mean diurnal variation 7-day moving windows (Falge et al., 2001). A linear regression model was fitted between daytime GPP and APAR for each 7-day moving window to estimate LUE, where LUE is the slope of the line.

Hyperspectral HCRF measurements and NDSI estimates
Ground surface HCRF spectra were measured every 15 min between sunrise and sunset from 15 July 2011 until 31 December 2012 using two FieldSpec 3 spectrometers with fibre optic cables (Table 1) (ASD Inc., Colorado, USA). The spectroradiometers cover the spectral range from 350 to 1800 nm and have a FOV of 25 • . The spectral resolution is 3 nm at 350-1000 nm and 10 nm at 1000-1800 nm and the sampling interval is 1.4 nm at 350-1000 nm and 2 nm at 1000-1800 nm. From these data, 1 nm spectra were calculated by using cubic spline interpolation functions. One sensor head was mounted on a rotating head 10.5 m above the surface (at the same tower including instruments to measure meteorological variables), providing measurements of the herbaceous vegetation from seven different viewing angles in a transect underneath the tower (nadir, 15, 30, and 45 • off-nadir angles towards east and west). No trees or effects of shading of trees are present in the IFOV of the data used in this study (Fig. 1).
A reflective cosine receptor is used to measure full-sky irradiance by having the second sensor head mounted on a 2 m high stand pointing to a Spectralon panel (Labsphere Inc., New Hampshire, USA) under a glass dome. Each sensor measurement starts with an optimisation to adjust the sensitivity of the detectors according to the specific illumination conditions at the time of measurement. The optimisation is followed by a dark-current measurement to account for the noise generated by the thermal electrons within the ASD instruments that flow even when no photons are entering the device. The measurement sequence starts with a full-sky-irradiance measurement, followed by measurements of the seven angles of the land surface and finalised by a second full-sky-irradiance measurement. Thirty scans are averaged to one measurement to improve the signal-to-noise ratio for each measurement (optimisation, dark current, fullsky irradiance, and each of the seven target measurements). The full measurement sequence takes less than 1 min. The two ASD instruments are calibrated against each other before and after each rainy season. Poor-quality measurements caused by unfavourable weather conditions, changing illumination conditions, and irregular technical issues were filtered by comparing full-sky solar irradiance before and after the target measurements (Huber et al., 2014). The spectral HCRF was derived by estimating the ratio between the ground surface radiance and full-sky irradiance. For a complete description/illustration of the spectroradiometer setup, the measurement sequence, and the quality control, see Huber et al. (2014).
where ρ i and ρ j are the daily median HCRF in two separate single wavelengths ( i and j ) between 350 and 1800 nm. In order to minimise the influence of errors we used daily median hyperspectral HCRF in the analysis (since the median provides the most common model output and is thereby more robust against outliers than average values).

Effects of varying sun and sensor viewing geometry on NDSI
The effects of variable solar zenith angles on different NDSI combinations were studied with nadir HCRF measurements. In order to capture the seasonal dynamics, data were taken over 15 days during four periods: (1)  Only days with full data coverage were used in order not to include bias in the results from days with incomplete data sets. The median HCRF of the 15 days was calculated for each wavelength for every 15 min between 08:00 and 18:00 (UTC). These HCRF values were combined into NDSI with different wavelength combinations. Finally, daily mean and standard deviation for all wavelength combinations were calculated. Diurnal variability in the NDSI was assessed with the coefficient of variation (COV), which is the ratio between the standard deviation and the mean. The COV gives an indication of effects related to variable solar zenith angles.
To capture directional effects in the NDSI related to the variable view zenith angles (15, 30, and 45 • off-nadir angles towards east and west) the NDSI was calculated using median HCRF values from the four above-mentioned periods for the different viewing angles. Only data measured between 12:00 and 14:00 (UTC) were used so as to avoid effects of variable solar zenith angles. The anisotropy factor (ANIF) is defined as the fraction of a reflected property at a specific view direction relative to the nadir, and it was calculated by where NDSI(λ, θ ) is NDSI for the different wavelengths (λ) and the different viewing angles (θ ), and NDSI 0 (λ) is the nadir-measured NDSI (Sandmeier et al., 1998).

Relationship between hyperspectral HCRF, NDSI, and ecosystem properties
We examined the relationship between predictor variables (daily median hyperspectral HCRF, and NDSI from nadir observations) and response variables (biomass, GPP, LUE, and FAPAR). A comparison between fitted linear and exponential regression models indicated no improvement by fitting exponential regression models; we hence choose to use linear regression analysis (Supplement). Possible errors (random sampling errors, aerosols, dust or water on the sensor heads, electrical senor noise, filtering and gapfilling errors, errors in correction factors, sensor drift, and instrumentation errors) can be present in predictor and response variables. We thereby used a reduced major axis linear regression to account for errors in both the predictor and response variables when fitting the regression lines. In order to estimate the robustness of the empirical relationships, we used a bootstrap simulation methodology, where the data sets were copied 200 times (Richter et al., 2012). The runs generated 200 sets of slopes, intercepts, and coefficients of determination (R 2 ), from which the median and standard deviation were estimated. The generated statistical models were validated against the left-out subsamples within the bootstrap simulation method by calculating the root-mean-square error (RMSE) and the relative RMSE (RRMSE = 100 × RMSE × mean(observed) −1 ); median and standard deviation were estimated. Within the regression analysis, all variables used were repeated observations of the same measurement plot. The dependent and independent variables are hence temporally autocorrelated and cannot be regarded as statistically independent. We thereby choose not to present any statistical significance. The analyses, however, still indicate how closely coupled the explanatory variables are with the ecosystem properties. A filter was created for the analysis between NDSI and ecosystem properties; all NDSI combinations with a COV higher than 0.066 in any of the four periods (dry season, fast growth period, peak of the growing season, and senescent period) and all NDSI combinations with ANIF values higher than 1.2 and lower than 0.8 in any of the four periods were filtered. The ANIF thresholds of 1.2 and 0.8 and the COV threshold of 0.066 were used since values then vary less than 20 % due to effects of variable sun-sensor geometry. NDSI including the water absorption band (1300-1500 nm) was also removed as it is strongly sensitive to atmospheric water content and is less suitable for spatial extrapolation of ecosystem properties using air/spaceborne sensors (Asner, 1998). Finally, NDSI combinations including wavelengths between 350 and 390 nm were removed owing to low signalto-noise ratio in the ASD sensors (Thenkabail et al., 2004).

Seasonal dynamics in meteorological variables, ecosystem properties, and hyperspectral HCRF
Daily average air temperature at 2 m height ranged between 18.4 and 37.8 • C, with low values during winter and peak val- ues at the end of the dry season (Fig. 2a). Yearly rainfall was 486 and 606 mm for 2011 and 2012, respectively. Soil moisture ranged between 1.9 and 14.1 %, and it clearly followed the rainfall patterns ( Fig. 2b and c). The CO 2 fluxes were low during the dry period and high during the rainy season (July-October) (Fig. 2e). The LUE followed GPP closely (Fig. 2f). FAPAR was low at the start of the rainy season, followed by a maximum towards the end of the rainy season, and then slowly decreased over the dry season (Fig. 2g). The range in HCRF is large across the spectral space, and would hide the seasonal dynamics in hyperspectral HCRF if directly shown. Therefore, to clearly illustrate these seasonal dynamics, the ratio between daily median nadir HCRF and the average HCRF for the entire measurement period was calculated for each wavelength (350-1800 nm). This gives a fraction of how the HCRF for each wavelength varies over Figure 3. The coefficient of variation (COV), i.e. the ratio between daily standard deviation and the daily mean (measurements taken between 08:00 and 18:00 (UTC)), for different normalised difference spectral index (NDSI) wavelength ( i, j ) combinations for 12 days at the peak of the growing season 2011 (day of year 237-251; n = 576). The COV indicates how strongly the NDSI are affected by variable sun angles. The upper right half of the chart shows the unfiltered R 2 values, whereas the lower left half shows filtered R 2 , based on the filtering criteria described in Sect. 2.6. the measurement period in relation to the average of the entire period (Fig. 2d). In the visible (VIS) part of the spectrum (350-700 nm) there was a stronger absorption during the second half of the rainy season and at the beginning of the dry season than during the main part of the dry season and the start of the rainy season. There was stronger NIR absorption (700-1300 nm) at the end of the rainy season and the beginning of the dry season, whereas the absorption decreased along with the dry season. Strong seasonal variation was observed in the water absorption region around 1400 nm following the succession of rainy and dry seasons. HCRF in the SWIR (1400-1800 nm) generally followed the seasonal dynamics of the visible part of the spectrum.

Effects of sensor viewing geometry and variable sun angles on NDSI
The strongest effects of solar zenith angles and variable viewing geometry on NDSI were observed at the peak of the growing season 2011 (Figs. 3 and 4, and S1-S5 in the Supplement). In the main section of the paper, we hence choose to present the figures from this period; figures from remaining periods are located in the Supplement. The most pronounced effects of solar zenith angles were observed for NDSI combining SWIR and NIR wavelengths, and with VIS wavelengths between 550 and 700 nm (n = 576) (Fig. 3). The same effects were seen for the view zenith angles; the strongest effects were seen for NDSI with SWIR and NIR combinations, as well as for VIS wavelengths between 550 and 700 nm (Fig. 4). Remaining VIS wavelengths were less affected. It was also clear that ground surface anisotropy increased strongly as a function of increasing viewing angle (Fig. 4). Moreover, some band combinations showed already angular sensitivity at view zenith angles of 15 • , while other band combinations only manifest anisotropic behaviour with higher view angles. Some band combinations, however, do not show any increased anisotropy at all (areas coloured in green in all six plots).

Biomass
HCRF values for all wavelengths except the water absorption band at 1100 nm were strongly correlated with biomass . The water absorption band (1300-1500 nm) was removed as it is strongly sensitive to atmospheric water content, and wavelengths between 350 and 390 nm were removed owing to low signal-to-noise ratio in the ASD sensors. (Fig. 5a). The strongest correlation was found at ρ 1675 (median ±1 SD; r = −0.88 ± 0.09), but biomass was almost equally well correlated with blue, red, and NIR wavelengths. All presented correlations and relationships throughout the text are based on filtered data. Negative correlations indicate that the more biomass there is, the higher the absorption, and hence the lower the HCRF. A small peak of pos- itive correlation is seen at 1120-1150 nm. NDSI combinations with HCRF in the red edge (ρ 680 -ρ 750 ) and HCRF in the VIS region explained seasonal dynamics in biomass well (Fig. 6a). The strongest relationship (R 2 = 0.88 ± 0.07; RRMSE = 18.6 ± 5.7 %) between NDSI and biomass was found for NDSI combining 705 and 587 nm (NDSI[705, 587]) (Table 2, Fig. 7a).

Gross primary productivity
The relationship between GPP and nadir-measured hyperspectral HCRF is inverted as compared to other correlation coefficient lines (Fig. 5b), since GPP is defined as a withdrawal of CO 2 from the atmosphere with higher negative values for a larger CO 2 uptake. The seasonal dynamics in GPP was strongly positively correlated with HCRF in the blue, red, SWIR wavelengths, and the water absorption band at 1100 nm, whereas it was strongly negatively correlated with the NIR HCRF. The study revealed the strongest positive and negative correlations for HCRF at 682 nm (r = 0.70 ± 0.02) and 761 nm (r = −0.74 ± 0.02), respectively. NDSI combinations that explained most of the GPP variability were different combinations of the VIS and NIR or red and SWIR wavelengths (Fig. 6b). However, the strongest relationship was seen at NDSI[518, 556] (R 2 = 0.86 ± 0.02; RRMSE = 34.9 ± 2.3 %) (Table 2; Fig. 7b).

Light use efficiency
LUE was negatively correlated with HCRF in the blue and red spectral ranges and in the water absorption band at 1100 nm, and it was positively correlated in the NIR wavelengths (Fig. 5c). HCRF at 761 nm yielded the strongest positive correlation (r = 0.87 ± 0.01). When combining the different wavelengths to NDSI, the VIS wavelengths explained variation in LUE well, with the strongest relationships in the red and blue parts of the spectrum (Fig. 6c). LUE correlated most strongly with NDSI[436, 688] (R 2 = 0.81 ± 0.02; RRMSE = 52.8 ± 3.8 %)) ( Table 2; Fig. 7c).

Fraction of photosynthetically active radiation absorbed by the vegetation
FAPAR was negatively correlated with nadir-measured HCRF for most wavelengths (Fig. 5d); the higher the FA-PAR, the higher the absorption, and thereby the lower the HCRF. The strongest correlation was found at a blue wavelength ρ 412 (r = −0.9 2± 0.01). When wavelengths were combined to NDSI, combining violet/blue with NIR and SWIR wavelengths generated the NDSI with the strongest relationships (Fig. 6d) with a maximum R 2 of 0.81 ± 0.02 (RRMSE = 14.6 ± 0.7 %) for NDSI [399,1295] (Table 2; Fig. 7d).

4630
T. Tagesson et al.: Deriving seasonal dynamics in ecosystem properties of semi-arid savanna grasslands Table 2. Wavelengths of the hemispherical conical reflectance factors (HCRF) (ρ i, j ; nm) used in the normalised difference spectral indices (NDSI) that generated the strongest correlations with ecosystem properties. DW is dry weight, FAPAR is the fraction of photosynthetically active radiation absorbed by the vegetation, Avg. is average, SD is standard deviation, RMSE is root-mean-square error, and RRMSE is relative RMSE.

Effects of sensor viewing geometry and variable sun angles on the NDSI
Effects of solar zenith angles and sensor viewing geometry were similar (Figs. 3 and 4), since they affect HCRF measurements in a similar way (Kimes, 1983). In dense and erectophile canopies, HCRF increases with sensor viewing and solar zenith angles, because a larger fraction of the upper vegetation canopy is viewed/illuminated, whereas the shadowed lower part of the canopy contributes less to the measured signal as shown previously by several studies (Huete et al., 1992(Huete et al., , 2014Jin et al., 2002;Kimes, 1983). However, the radiative transfer within a green canopy is complex, and differs across the spectral region (Huber et al., 2014). Less radiation is available for scattering of high-absorbance spectral ranges (such as the VIS wavelengths), and this tends to increase the contrast between shadowed and illuminated areas for these wavelengths, whereas in the NIR and SWIR ranges, more radiation is scattered and transmitted, which thereby decreases the difference between shadowed and illuminated areas within the canopy (Kimes, 1983;Hapke et al., 1996). A recognised advantage of NDSI calculations is that errors/biases that are similar in both wavelengths included in the index are suppressed by the normalisation. However, for a given situation where errors/biases are different for the wavelengths used, such as effects generated by sun-sensor geometry, errors/biases will affect the value of the index. This was also the case at the Dahra field site: NDSI values were strongly affected at wavelength combinations with large differences in effects of variable solar zenith angles ( Fig. 6 in Huber et al., 2014) and variable view zenith angles ( Fig. 4 in Tagesson et al., 2015b). This effect is especially pronounced in the case of low index values (closer to 0), whereas larger index values (closer to 1 and −1) become less sensitive. The relative HCRF difference between NIR and SWIR is lower compared to indices including the VIS domain; NIR/SWIR-based indices thereby generate lower NDSI values with higher sensitivity to sun-sensor geometry-generated differences between included wavelengths (Figs. 3 and 4). This can also be seen in the SIWSI-NDVI comparison by Huber et al. (2014); SIWSI had large relative differences depending on viewing angle (∼ 70 %), whereas NDVI had relatively small differences (∼ 5 %) ( Fig. 10 in Huber et al., 2014). Fensholt et al. (2010a) showed the same to be true in a comparison between SIWSI and NDVI based on MODIS data: SIWSI was insensitive to day-to-day variations in canopy water status due to effects of solar zenith angles and sensor viewing geometry blurring the signal. A strong diurnal dynamic does not necessarily mean a poor NDSI. For example, the photochemical reflectance index (PRI) was created for assessing diurnal dynamics in the xanthophyll cycle activity (Gamon et al., 1992). Stomatal closure due to high temperatures could also influence diurnal dynamics of vegetation properties (Lasslop et al., 2010), and hence the diurnal dynamics of NDSI. However, diurnal variation in reflectance caused by diurnal variability in vegetation status is assumed minor in relation to the diurnal variability caused by changes in solar zenith angles. Additionally, in our study we are interested in relationships in seasonal dynamics between ecosystem properties and NDSI; diurnal variation can thereby interfere and introduce uncertainty into such relationships.
The importance of directional effects for the applicability of normalised difference spectral indices has been pointed out as an issue in numerous papers (e.g. Holben and Fraser, 1984;van Leeuwen et al., 1999;Cihlar et al., 1994;Fensholt et al., 2010b;Gao et al., 2002). This study confirms these challenges for NIR/SWIR-based indices, but the results also indicate several wavelength combinations from which these effects are less severe and potentially applicable to EO data without disturbance from viewing/illumination geometry for this type of vegetation. Multi-angular HCRF data provide additional information of, for example, canopy structure, photosynthetic efficiency, and capacity (Bicheron and Leroy, 2000;Asner, 1998;Pisek et al., 2013;Huber et al., 2010), and this unique in situ-based, multi-angular, hightemporal-resolution data set may thus be used for future research of canopy radiative transfer and BRDF (bidirectional reflectance distribution function) modelling (Jacquemoud et al., 2009;Bicheron and Leroy, 2000). The multi-angular data set is also highly valuable for evaluation and validation of satellite-based products, where the separation of view angle and atmospheric effects can only be done using radiative transfer models (Holben and Fraser, 1984).

Biomass
The strong correlation between biomass and most of the spectrum indicates the strong effects of phenology on the seasonal dynamics in the HCRF spectra (Fig. 5a). Variability in VIS (350-700 nm) HCRF for vegetated areas is strongly related to changes in leaf pigments (Asner, 1998), and this can also be seen in Fig. 2d since absorption was much stronger during the rainy (growing) season than during the dry season. Previous studies have generally shown positive relationships between NIR HCRF and biomass since a large fraction of NIR radiation is reflected in green healthy vegetation to avoid overheating (e.g. Hansen and Schjoerring, 2003;Asner, 1998). Here, a strong negative relationship between NIR HCRF and dry weight biomass is generally observed (Fig. 5a), indicating stronger NIR absorption with increased biomass. However, a strong positive NIR HCRF correlation with vegetation water content was seen (figure not shown). A possible explanation could be that the sampled biomass at the end of the rainy season contained some senescent vegetation, and a correlation against vegetation water content is hence closer to green healthy vegetation. This relationship is, however, interesting and should be studied further to better understand the respective importance of canopy water and leaf internal cellular structure for the NIR HCRF of herbaceous vegetation characterised by erectophile leaf angle distribution in semi-arid regions. We found the strongest correlation for biomass with a SWIR wavelength, thereby confirming the studies by  and Psomas et al. (2011) in that SWIR wavelengths are good predictors of LAI or biomass. The NDVI is known to saturate at high biomass because the absorption of red light at ∼ 680 nm saturates at higher biomass loads, whereas the NIR HCRF continues to increase due to multiple scattering effects (Mutanga and Skidmore, 2004;Jin and Eklundh, 2014). Several studies have shown that NDSI computed with narrowband HCRF improve this relationship by choosing a wavelength region not as close to the maximum red absorption at ∼ 680 nm, for example using shorter and longer wavelengths of the red edge (700-780nm) (Cho et al., 2007;Mutanga and Skidmore, 2004; and NIR and SWIR wavelengths (Psomas et al., 2011;. The NDSI with the strongest correlation to biomass was computed using red-edge HCRF (ρ 705 ) and green HCRF (ρ 587 ). Vegetation stress and information about chlorophyll and nitrogen status of plants can be extracted from the red-edge region (Gitelson et al., 1996). Wavelengths around ρ 550 are located right at the peak of green reflection and closely related to the total chlorophyll content, leaf nitrogen content, and chlorophyll / carotenoid ratio and have previously been shown to be closely related to biomass (Inoue et al., 2008;Thenkabail et al., 2012).

Gross primary productivity
The maximum absorption in the red wavelengths generally occurs at 682 nm as this is the peak absorption for chlorophyll a and b (Thenkabail et al., 2000), and this was also the wavelength that was most strongly correlated with GPP. HCRF at 682 nm was previously shown to be strongly related to LAI, biomass, plant height, NPP, and crop type discrimination (Thenkabail et al., 2004. The NDSI with the strongest relationship to GPP was based on HCRF in the vicinity of the green peak. The PRI normalises HCRF at 531 and 570 nm and it was suggested for detection of diurnal variation in the xanthophyll cycle activity (Gamon et al., 1992), and it is commonly used for estimating productivity efficiency of the vegetation (e.g. Soudani et al., 2014). The present study thereby confirms the strong applicability of the wavelengths in the vicinity of the green peak for vegetation productivity studies. Again, wavelengths around the green peak are related to the total chlorophyll content, leaf nitrogen content, chlorophyll / carotenoid ratio, and biomass (Inoue et al., 2008;Thenkabail et al., 2012).

Light use efficiency
Both LUE and GPP were most strongly correlated with HCRF at 761 nm, which is the oxygen A band within the NIR wavelengths. HCRF at 761 nm is commonly used for estimating solar-induced chlorophyll fluorescence due to radiation emitted by the chlorophyll, and it has been suggested as a direct measure of health status of the vegetation (Meroni et al., 2009). Earth observation data for estimating fluorescence should have very high spectral resolution (< 10 nm), but considering the rapid technical development within sensors for hyperspectral measurements, fluorescence possibly has strong practical potential for monitoring vegetation status (Meroni et al., 2009;Entcheva Campbell et al., 2008). Globally mapped terrestrial chlorophyll fluorescence retrievals are already produced from the GOME-2 instrument at a spatial resolution of 0.5 • × 0.5 • , but hopefully this will also be available with EO sensors of higher spatial and temporal resolution in the future (Joiner et al., 2013).
The strongest wavelength combinations for estimating LUE for this semi-arid ecosystem was NDSI [688,435]. The 688 nm wavelength is just at the base of the red-edge region, again indicating the importance of this spectral region for estimating photosynthetic activity. The wavelength at 435 nm is at the centre of the blue range characterised by chlorophyll utilisation, and strongly related to chloro-4632 T. Tagesson et al.: Deriving seasonal dynamics in ecosystem properties of semi-arid savanna grasslands phyll a and b, senescing, carotenoid, loss of chlorophyll, and vegetation browning (Thenkabail et al., 2004. The NDSI [688,435] thereby explores the difference between information about chlorophyll content and information about senescence of the canopy, which should be a good predictor of ecosystem-level photosynthetic efficiency.

Fraction of photosynthetically active radiation absorbed by the vegetation
FAPAR is an estimate of radiation absorption in the photosynthetically active spectrum and thereby strongly negatively correlated with most parts of the spectrum (Fig. 5d). FA-PAR remained high during the dry season because of standing dry biomass that was slowly degrading over the dry season (Fig. 2g). The seasonal dynamics in FAPAR is thereby strongly related to senescence of the vegetation, which explains why FAPAR was most strongly correlated with blue wavelengths (ρ 412 ). Several studies have reported a strong relationship between NDVI and FAPAR (e.g. Tagesson et al., 2012;Myneni and Williams, 1994;Fensholt et al., 2004), but this relationship has been shown to vary for the vegetative phase and the periods of senescence (Inoue et al., 1998;Tagesson et al., 2015b). As showed by Inoue et al. (2008), and confirmed by this study, new indices combining blue with NIR wavelengths could be used for estimating FAPAR for the entire phenological cycle. This result has implications for studies using the LUE approach for estimating C assimilations (Hilker et al., 2008).

Outlook and perspectives
Very limited multi-angular hyperspectral in situ data exist, even though they have been, and will continue to be, extremely valuable for an improved understanding of the interaction between ground surface properties and radiative transfer. In this study, we have presented a unique in situ data set of multi-angular, high temporal resolution hyperspectral HCRF (350-1800 nm) and demonstrated the applicability of hyperspectral data for estimating ground surface properties of semi-arid savanna ecosystems using NDSI. The study was conducted in spatially homogeneous savanna grassland, suggesting that the results should be commonly applicable for this biome type. However, attention should be paid to sitespecific details that could affect the indices, such as species composition, soil type, biotic and abiotic stresses, and stand structure. Additionally, the biophysical mechanisms behind the NDSIs are not well understood at the moment, and further studies are needed to examine the applicability of these indices to larger regions and other ecosystems. Due to it being a two-band ratio approach, NDSI does not take full advantage of exploring the rich information given by the hyperspectral HCRF measurements. In the future, this hyperspectral HCRF data set could be fully explored using, for example, derivative techniques; multivariate methods; and creation, parameteri-sation, and evaluation of BRDF and radiative transfer models. Even though several other methods exists which fully exploit the information in the hyperspectral spectrum, results of the present study still indicate the strength of normalised difference indices for extrapolating seasonal dynamics in properties of savanna ecosystems. A number of wavelength spectra that are highly correlated with seasonal dynamics in properties of semi-arid savanna ecosystems have been identified. The relationships between NDSI and ecosystem properties were better determined than, or at the same level as, results of previous studies exploring relationships between hyperspectral reflectance and ecosystem properties (Kumar, 2007;Cho et al., 2007;Mutanga and Skidmore, 2004;Psomas et al., 2011;Ide et al., 2010). By also studying the impact from varying viewing and illumination geometry, the feasibility and applicability of using indices for upscaling to EO data were evaluated. As such, the results presented here offer insights for assessment of ecosystem properties using EO data, and this information could be used for designing future sensors for observation of ecosystem properties of the Earth.
The Supplement related to this article is available online at doi:10.5194/bg-12-4621-2015-supplement.