Monitoring vegetation condition using microwave remote sensing: the standardized vegetation optical depth index (SVODI)

Vegetation conditions can be monitored on a global scale using remote sensing observations in various wavelength domains. In the microwave domain, data from various spaceborne microwave missions are available from the late 1970s onwards. From these observations, vegetation optical depth (VOD) can be estimated, which is an indicator of the total canopy water content and hence of above-ground biomass and its moisture state. Observations of VOD anomalies would thus complement indicators based on visible and near-infrared observations, which are primarily an indicator of an ecosystem's photosynthetic activity.

Reliable long-term vegetation state monitoring needs to account for the varying number of available observations over time caused by changes in the satellite constellation. To overcome this, we introduce the standardized vegetation optical depth index (SVODI), which is created by combining VOD estimates from multiple passive microwave sensors and frequencies. Different frequencies are sensitive to different parts of the vegetation canopy. Thus, combining them into a single index makes this index sensitive to deviations in any of the vegetation parts represented. SSM/I-, TMI-, AMSR-E-, WindSat- and AMSR2-derived C-, X- and Ku-band VODs are merged in a probabilistic manner resulting in a vegetation condition index spanning from 1987 to the present.

SVODI shows similar temporal patterns to the well-established optical vegetation health index (VHI) derived from optical and thermal data. In regions where water availability is the main control on vegetation growth, SVODI also shows similar temporal patterns to the meteorological drought index scPDSI (self-calibrating Palmer drought severity index) and soil moisture anomalies from ERA5-Land. Temporal SVODI patterns relate to the climate oscillation indices SOI (Southern Oscillation index) and DMI (dipole mode index) in the relevant regions. It is further shown that anomalies occur in VHI and soil moisture anomalies before they occur in SVODI.

The results demonstrate the potential of VOD to monitor the vegetation condition, supplementing existing optical indices. It comes with the advantages and disadvantages inherent to passive microwave remote sensing, such as being less susceptible to cloud coverage and solar illumination but at the cost of a lower spatial resolution.

The index generation is not specific to VOD and could therefore find applications in other fields.

The SVODI products (Moesinger et al., 2022) are open-access under Attribution 4.0 International and available at Zenodo, https://doi.org/10.5281/zenodo.7114654.

1 Introduction

Monitoring vegetation conditions by remote sensing is important for a variety of purposes, such as agricultural yield prediction (Petersen, 2018; Crocetti et al., 2020), forestry (Pause et al., 2016) and fire ecology (Szpakowski and Jensen, 2019) and to track long-term ecosystem changes (Vogelmann et al., 2012). On a global scale, observing vegetation conditions from space allows for cost-effective, long-term and spatially consistent analyses.

Numerous variables and metrics have been developed to monitor vegetation conditions from spaceborne observations. Some use the spectral or radiometric information directly to create a feature related to the vegetation. This includes features such as the normalized difference vegetation index (NDVI), which is widely used as a measure of live green vegetation (Huang et al., 2021; Tucker et al., 2005), or the cross-polarization ratio (CR), which is related to polarization changes of active microwaves caused by vegetation structure and moisture content (Vreugdenhil et al., 2020).

Often, radiometric information is translated into biogeophysical or chemical variables such as the leaf area index (LAI) or the fraction of absorbed photosynthetic radiation (Dorigo et al., 2007). An ecosystem's condition can also be observed through its response to stress, e.g. by measuring changes in land surface temperature (Kogan, 1990) or evaporation (Martens et al., 2017).

All these features have in common that they show some aspect of the vegetation at a given time and location.

To assess whether the state of the vegetation is unusual at a given time and location, it is usually compared against the expected value at that time of year, derived from long-term observations. There are multiple ways to do this. The most straightforward way is to calculate anomalies by subtracting the multiyear seasonal average from the observation. This allows one to directly see whether an observation is higher or lower than usual. The drawback of such raw anomalies is that their magnitude depends on the average conditions at a given location. Therefore, the anomalies between different locations cannot be compared against each other and it requires expert knowledge to know whether an anomaly of a certain magnitude is a very strong outlier or a minor deviation (Katz and Glantz, 1986).

This can be solved by expressing the vegetation condition as an index. While anomalies show deviations as absolute differences to some mean, indices show the likelihood of observing a deviation of a certain magnitude. Indices are easier to interpret as they follow a well-defined distribution which allows one to discern quickly whether a value is relatively high or low. Some well-known example indices are the vegetation condition index (VCI, computed from NDVI; Kogan, 1990, 1997, 2001), the temperature condition index (TCI, computed from observations in the thermal domain; Kogan, 1990, 1997, 2001), and the standardized precipitation index (SPI, from precipitation estimates; McKee et al., 1993).

Over the past 4 decades, various platforms carrying multi-frequency microwave radiometers have been orbiting the Earth. From these observations it is possible to derive the vegetation optical depth (VOD), which describes the attenuation of microwave radiation by vegetation (Jackson and Schmugge, 1991; Meesters et al., 2005). The higher the vegetation water content and the shorter the wavelength, the more the vegetation attenuates the radiation (Jackson and Schmugge, 1991; Owe et al., 2008). Each frequency band is sensitive to slightly different parts of the vegetation, with short wavelengths such as those measured by the Ku-band being mostly related to the canopy top and leaves, while longer wavelengths are also sensitive to the woody vegetation parts (Owe et al., 2008; Rodríguez-Pérez et al., 2018). Compared to indices derived from optical data, VOD saturates less quickly for dense canopies and is therefore more sensitive to fluctuations in densely vegetated areas (Liu et al., 2015; Frappart et al., 2020). Among many other things, VOD has been used to analyse the vegetation's response to droughts in the Amazonian tropics (Liu et al., 2018) and in the Pannonian Basin (Crocetti et al., 2020) to determine deforestation in the tropics (van Marle et al., 2016) and estimate gross primary production (Teubner et al., 2019).

Long-term VOD datasets, such as VODCA (Moesinger et al., 2020) or the dataset presented by Liu et al. (2015), allow for monitoring vegetation conditions over decadal timescales. It might seem trivial to create an index from any of these datasets by calculating the seasonal anomalies and standardizing them. However, these datasets are based on averaging all available VOD values from different sensors. This causes the merged data to be heteroscedastic, where periods with fewer observations are noisier than periods with more observations due to the law of large numbers. The noise level of the sensors used also differs over time, where newer sensors generally are less noisy. High noise levels increase the probability of a value to be extreme, and therefore extreme values are more likely to occur in periods with few observations or noisy sensors. This hampers comparisons of extreme events over longer time periods.

To solve this issue, the standardized VOD index (SVODI) is proposed, which uses a probabilistic merging method to generate a long-term dataset for global vegetation condition monitoring based on VOD. After a technical evaluation, its relationship to other vegetation-related indices is explored. This assures that SVODI behaves reasonably in the case of an event affecting the vegetation and gives insight into how it differs from the currently used indices for vegetation condition monitoring.

2 Data

2.1 Vegetation optical depth datasets

2.1.1 The land parameter retrieval model (LPRM)

VOD estimates from various microwave radiometers and frequencies have been obtained with LPRM v6.0 (van der Schalie et al., 2017; Owe et al., 2008; Meesters et al., 2005), which is a forward radiative transfer model based on the work of Mo et al. (1982). It simulates the top-of-atmosphere brightness temperature for a wide range of surface conditions. It retrieves soil moisture and VOD analytically using polarized microwave data and Ka-band surface temperature estimations (Holmes et al., 2009) without relying on external information on the vegetation. LPRM assumes a thermal equilibrium between the surface and vegetation (Owe et al., 2008). Due to uneven solar heating during the day, daytime LPRM retrievals are still very experimental with an unknown error magnitude. In line with other studies using LPRM data (Dorigo et al., 2017; Moesinger et al., 2020), only nighttime VOD observations are used.

2.1.2 Sensor specifications

For this study, the exact same data as for VODCA (Moesinger et al., 2020) are used, namely VOD data derived from the radiometers SSM/I, TMI, AMSR-E, WindSat and AMSR2. An overview of the specifications can be found in Table 1. All sensors but TMI have a sun-synchronous circular orbit, resulting in global coverage.

SSM/I (Special Sensor Microwave/Imager) F08, F11 and F13 on board DMSP satellites are used for a total time span of 1987 to 2015. VOD retrieved from the Ku-band with a resolution of 69×43 km is used (Wentz, 1997).

TMI, the TRMM Microwave Imager on board TRMM, is used from 1997 to 2015. The VOD retrieved from the X- and Ku-band are used. Flying in a non-near-polar orbit, it only covered the area 35^∘ N–35^∘ S until 2001, when a boost in altitude increased it to about 37^∘ N–37^∘ S. This also decreased the spatial resolution of the X-band from 63×37 to 72×43 km and of the Ku-band from 30×18 to 35×21 km (Kummerow et al., 1998).

AMSR-E, the Advanced Microwave Scanning Radiometer for EOS (the Earth Observing System) on board Aqua, is used from 2002 to 2011. The VODs retrieved from the C-, X- and Ku-band are used, which have a spatial footprint of 75×43, 51×29 and 27×16 km, respectively. Only the descending overpass is used, which passes the Equator at 01:30 UTC (Knowles et al., 2006; Kawanishi et al., 2003).

WindSat on board Coriolis observes the C-, X- and Ku-band with a spatial resolution of 39×71, 25×38 and 16×27 km, respectively, and the retrieved VOD values from 2003 to 2012 are used (Gaiser et al., 2004). WindSat is still in orbit and functional, but no access to data past 2012 has been given.

AMSR2, the Advanced Microwave Scanning Radiometer 2 on board GCOM-W1, is used from June 2012 onwards. It is the follow-up to AMSR-E and as such is very similar but with slightly higher spatial resolution of 62×35, 42×24 and 22×14 km for the C-, X- and Ku-band, respectively. Another improvement is the addition of a second C-band (7.3 GHz) that can be used in case radio-frequency interference (RFI) affects the primary C-band (6.9 GHz) (Markus et al., 2018). For the C- and X-band, all retrieved daytime VODs are used. For the Ku-band, preliminary analysis revealed that the VOD retrievals after 1 August 2017 abruptly dropped globally, possibly due to a calibration issue. While the exact reason for the change is not known, the data are deemed unreliable and are not used after that date.

2.2 Auxiliary data

Multiple auxiliary datasets are used to evaluate SVODI. Most of these datasets do not follow a standard normal distribution, either by design or because the preprocessing (regridding and temporal resampling) changed their distribution. To facilitate the comparison of all datasets, they are standardized using a basic Z-score normalization:

$$\begin{array}{}\text{(1)}& x_{\mathrm{standardized}}={\displaystyle \frac{\left(x-{\mathit{\mu }}_x\right)}{{\mathit{\sigma }}_x}},\end{array}$$

where x_standardized is the standardized data; x is the original data; and μ_x and σ_x are its mean and standard deviation, respectively. This ensures that all datasets for the comparison have a mean of 0 and a standard deviation of 1.

2.2.1 Vegetation health index, vegetation condition index and temperature condition index

The vegetation health index (VHI) (Kogan, 1990, 1997, 2001) is derived from optical and thermal observations. The general concept is to combine water and temperature stress indices into a combined vegetation health index (Kogan, 2001). As such, VHI is a weighted average of the vegetation condition index (VCI) and temperature condition index (TCI),

$$\begin{array}{}\text{(2)}& \mathrm{VHI}=\mathit{\alpha }\mathrm{VCI}+\left(\mathrm{1}-\mathit{\alpha }\right)\mathrm{TCI},\end{array}$$

where the weight α is traditionally set to 0.5 (Kogan, 1997). VCI is derived from NDVI and as such contains information about the greenness of the vegetation:

$$\begin{array}{}\text{(3)}& \mathrm{VCI}=\mathrm{100}\cdot {\displaystyle \frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}}\end{array}$$

where NDVI, NDVI_min and NDVI_max are the smoothed weekly NDVI and its multiyear minimum and maximum NDVI, respectively. VCI has been used successfully for drought monitoring and assessing vegetation conditions (Kogan, 1997) and in VHI is assumed to account for water stress.

The TCI is defined similarly but based on land surface temperature (LST),

$$\begin{array}{}\text{(4)}& \mathrm{TCI}=\mathrm{100}\cdot {\displaystyle \frac{{\mathrm{LST}}_{\max }-\mathrm{LST}}{{\mathrm{LST}}_{\max }-{\mathrm{LST}}_{\min }}},\end{array}$$

where LST, LST_min and LST_max are the smoothed weekly land surface temperature, its multiyear minimum and its multiyear maximum, respectively. The TCI increased the accuracy of drought monitoring by accounting for temperature stress and has been used to analyse the role of temperature in droughts (Kogan, 1997). By combining it with VCI, both water and temperature stress are accounted for.

VHI, VCI and TCI derived from AVHRR NDVI (Tucker et al., 2005) and the thermal channel 4 are used. The data are available at https://www.star.nesdis.noaa.gov (last access: 30 October 2022). They are downsampled from the original 4 km resolution to match our 0.25^∘ grid and are standardized using Eq. (1). The data start in 1981, but only values after 1987, the start date of SVODI, are used in this study.

An alternative to AVHRR data might be the newer MODIS data. However, the long-term availability of AVHRR allows for comparisons over the whole duration of SVODI, while MODIS is only available after 2000. Additionally, neither MODIS's higher spectral resolution (not relevant for VCI calculation) nor its higher spatial resolution (the AVHRR resolution is already much higher than our 0.25^∘ grid) is of any benefit in our application. Further AVHRR NDVI and MODIS NDVI correlate very strongly with each other (Bédard et al., 2006; Gallo et al., 2005; Zeng et al., 2013), reinforcing that the results (in the overlapping periods) would not change much by using MODIS instead of AVHRR. Overall therefore AVHRR is more suitable for this study.

2.2.2 Self-calibrating Palmer drought severity index

The self-calibrating Palmer drought severity index (scPDSI) (Wells et al., 2004; Van Der Schrier et al., 2013; Aldred et al., 2021) is a widely used index to track meteorological, agricultural and hydrological aspects of drought. It is used to analyse the relation between SVODI and meteorological droughts. One of the main drawbacks of the original PDSI (Palmer, 1965) was its lack of spatial comparability due to fixed weights and factors, which was remedied in scPDSI by adjusting them to the local climate. scPDSI models the soil moisture using a bucket model involving evaporation, recharge, runoff, loss and their complementary potential values. This gives a measure of how extreme the water conditions are at a certain time and place, which is useful for monitoring water stress. The scPDSI data are available at https://crudata.uea.ac.uk/cru/data/drought/ (last access: 30 October 2022).

2.2.3 ERA5-Land

ERA5-Land is a reanalysis of the global atmosphere, land surface and ocean waves since 1950 (Muñoz-Sabater et al., 2021; Hersbach et al., 2020). ERA5-Land models a plethora of land variables on a sub-daily temporal resolution. Of interest for our study is the modelled soil moisture at different depths, which is used assess the relation between SVODI and soil moisture at different depths (Sect. 4.2.2). ERA5-Land is available from the Climate Data Store at https://doi.org/10.24381/cds.e2161bac, last access: 30 October 2022.

2.2.4 SOI and DMI

The Southern Oscillation index (SOI) (Allan et al., 1991; Allan, 1998) and dipole mode index (DMI) (Saji et al., 1999; Saji and Yamagata, 2003) are useful to monitor large-scale climate oscillations in the tropics. SOI is derived from the sea-level atmospheric pressure difference between Tahiti and Darwin, Australia, while DMI is calculated from the difference in the sea surface temperature anomaly between the west and south-eastern tropical Indian Ocean. Among other things, both of them are linked to precipitation anomalies in Australia and north-eastern Africa, regions where the vegetation is water-limited and where they are therefore also expected to be linked to the vegetation condition (Hashimoto et al., 2019; Martens et al., 2018).

Both are available at NOAA, at https://psl.noaa.gov/gcos_wgsp/Timeseries/SOI (last access: 30 October 2022) and https://psl.noaa.gov/gcos_wgsp/Timeseries/DMI (last access: 30 October 2022), respectively.

3 Methods

3.1 SVODI calculation

SVODI is computed from C-, X- and Ku-band VOD from multiple sensors. It is assumed that all sensors and bands are equally fit as an indicator of the vegetation condition but show different aspects of it. Ku-band is mostly sensitive to surface canopy leaves, while longer wavelengths are also affected by the woody part (Owe et al., 2008; Rodríguez-Pérez et al., 2018). Still, the seasonal VOD anomalies of the three bands correlate strongly with each other (Fig. 1) as established in previous studies (Moesinger et al., 2020), which suggests that they behave similarly in case of a vegetation disturbance. To maximize the information contained in a microwave-based vegetation condition index, it makes sense to combine the information contained in all bands. Also, the high number of observations per day available due to the many sensors and bands is expected to yield an index robust to noise. Additionally, none of the bands span the whole time period (C-band VOD 2002 to present, X-band VOD 1997 to present, Ku-band VOD 1987 to 2017), so by merging them the longest possible time span is achieved.

Figure 1 Exemplary correlation coefficients between C-, X- and Ku-band VOD anomalies of AMSR-E, derived with LPRM (Sect. 2.1.1). Similar results are obtained for other sensors (not shown).

L-band VOD is notably absent from the list of used frequencies. It has very different temporal characteristics than the other bands due to it being mostly sensitive to slow structural changes in the vegetation (Konings et al., 2021), while SVODI is only concerned with short-term fluctuations. Therefore L-band VOD is not used.

3.1.1 Theory

Previous multivariate indices (Hao and AghaKouchak, 2013; Guo et al., 2019) have been constructed by first fitting a multivariate distribution with cumulative joint probability $P(X_{\mathrm{1}}\le x_{\mathrm{1}},\mathrm{\dots },X_n\le x_n)=p_{\mathrm{combined}}$ to n individual indices X_i, where $i\in \mathrm{1},\mathrm{\dots },n$. Then, the p_combined of each observation $(x_{\mathrm{1}},\mathrm{\dots },x_n)$ is transformed to the index by applying the standard normal percent point function (PPF, the inverse of the cumulative distribution function) to it.

But there is an issue with this approach: p_combined is not uniformly distributed between 0 and 1. Instead, it has a bias towards low values. This is illustrated for a theoretical bi-variate case with a Gaussian copula in Fig. 2: if some x₁, x₂ pairs are drawn and p is calculated for them, the distribution of p_combined is clearly not uniformly distributed. This causes the final index to have a negative bias. For example, consider the case $x_{\mathrm{1}}=x_{\mathrm{2}}=\mathrm{0}$, marked with red lines in Fig. 2. Instinctively, since both input indices show usual conditions, one would expect that the merged index should also show usual conditions. However, as seen in the figure, $P(X_{\mathrm{1}}\le \mathrm{0},X_{\mathrm{2}}\le \mathrm{0})<\mathrm{0.50}$, and therefore PPF(p_combined)<0. Therefore, even if all input indices are 0, the merged index will be negative. For a higher number of concurrent indices used as input, this effect is even stronger. In general, this would lead to an index with many more extreme negative events than positive ones (Fig. 2b).

Figure 2 Bi-variate example of a probabilistic multivariate index as in Hao and AghaKouchak (2013) and Guo et al. (2019) without scaling of p where $X_{\mathrm{1}},X_{\mathrm{2}}\sim N(\mathrm{0},\mathrm{1})$ and the covariance between them is 0.5. PDF (a); distribution of the CDF of samples from it (b); the resulting distribution of the index (c). The red lines mark an example case where $x_{\mathrm{1}}=x_{\mathrm{2}}=\mathrm{0}$ at different processing steps.

The negative bias of multivariate indices computed as above makes the index hard to interpret as the resulting index is no longer normally distributed. Also, in the case that the individual indices have data gaps, the expected mean depends on the available input indices, leading to a higher expected value for periods where fewer sensors are available than for periods with more sensors available. This issue is solved by scaling p_combined to a uniform distribution whose properties do not change depending on the number of available input indices. Therefore SVODI has no bias and can also be calculated if not all input datasets are available.

3.1.2 Implementation

After some basic preprocessing (temporal resampling of swath data to daily values and masking invalid values, the same as in Moesinger et al., 2020), the VOD values of each sensor and band, VOD_s,b, are transformed independently into standard normally distributed indices using the following workflow, which is independently applied for each band.

Long-term VOD changes are related to biomass changes (Frappart et al., 2020). The input datasets are therefore linearly detrended causing anomalies to correspond to deviations in the vegetation condition and not to long-term structural changes. Linear detrending might not be sufficient in areas experiencing rapid vegetation changes over a short time, such as deforestation, and climatologically unexpected biomass fluctuations, such as out-of-season harvests, as both of those cases would be registered by SVODI. However, a mitigating factor in both cases is the low spatial resolution of SVODI, which causes only very large-scale events to be registered. The main reason for using linear regression is that its simplicity guarantees that no deviations of interest are removed.

Table 1 Overview of VOD datasets used in this study with their temporal coverage, local ascending equatorial crossing times (AECT), whether the ascending (A) or descending (D) overpass is used, and frequencies (GHz) used for each product. The C- and X-band retrievals are based on van der Schalie et al. (2017) and the Ku-band retrievals on Owe et al. (2008).

Then, all VOD values of the respective band, VOD_b, are scaled to the values of the corresponding AMSR-E band to correct for bias between the different VOD_b values using improved piecewise linear cumulative distribution function (CDF) matching as described in Moesinger et al. (2020). AMSR-E is used as reference because of its high-quality observations, global availability and temporal overlap with most other sensors, a combination of features none of the other sensors possess. In detail, SSM/I, TMI and WindSat are matched to AMSR-E using temporally overlapping observations. AMSR2 does not have a temporal overlap with AMSR-E, and therefore direct matching using overlapping observations is not possible. Instead, between 37^∘ N and 37^∘ S AMSR2 is matched to TMI observations that were first scaled to AMSR-E. Beyond 37^∘ N and 37^∘ S, AMSR2 is directly matched to AMSR-E using the last 2 years of AMSR-E and first 2 years of AMSR2 to determine the scaling parameters.

The scaled VOD values are standardized in the following way: for a day of the year (DOY_i, where $i\in \mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{366}$), all VOD_b values from July 2002 to June 2017 that are less than 16 d from DOY_i are used to build an empirical distribution. This window size of 31 d was empirically chosen as a compromise between having enough values to build stable scaling parameters and the values not being biased in respect to the window centre due to the progressing seasonal VOD signal. The period July 2002 to June 2017 is chosen because all three frequency bands have observations during this period. Then, CDF scaling parameters are calculated to transform the empirical distribution to an N(0,1) distribution. Using these parameters, the VOD_b values at DOY_i are scaled. This is repeated for all DOYs and done independently for each band, resulting in individual indices X_s,b for each VOD_s,b.

Figure 3 Temporal subset of an example time series at different processing stages in Western Australia (24.9^∘ S, 125.625^∘ E), located in the Gibson Desert Nature Reserve. The vegetation is shrubland with sparse trees. (a) Original VOD time series. (b) The same VOD series after bias correction. Note that only the bias within each band is corrected and not to a common reference for all bands. (c) The indices created from each sensor and band as well as SVODI. (d) The number of observations contributing to a SVODI value for each day.

Figure 4 Fraction of pixels that are extreme of SVODI (a) and non-probabilistic merge of indices (b) over time.

The indices of the individual sensors and bands, X_s,b, are then joined by constructing a multivariate normal distribution with a zero mean. As discussed before, if $P(X_{\mathrm{1}}\le x_{\mathrm{1}},\mathrm{\dots },X_n\le x_n)=p_{\mathrm{combined}}$, then p_combined is not uniformly distributed but is biased towards low values. Therefore, p_combined is scaled to a uniform distribution. Technically this is done by drawing random samples from P and constructing an empirical CDF, CDF_P, from them. p_combined is then scaled to a uniform distribution by CDF_P(p_combined). The PPF is then applied to the scaled p_combined, resulting in a normally distributed index. The scaling of p_combined is numerically slightly unstable and can lead in very rare cases to extremely low SVODI values. Due to the computer's limited precision, these SVODI values always have the exact same value, i.e. −8.14. Therefore, observations where SVODI is lower than −8 are removed.

Since p_combined is always scaled to a uniform distribution, the number of available indices at a given time and location does not affect the distribution of SVODI and a continuous index starting in 1987 to the present is generated. For each SVODI value, a flag indicates which sensors and bands contributed to the final value.

3.2 Evaluation methods

SVODI describes the vegetation condition with regard to abnormal vegetation water content. There is no absolute reference to compare SVODI to, and therefore it is evaluated by comparing it to other well-established vegetation indicators. This is mostly done by basic correlation analysis but also by studying temporal shifts and the evolution of extreme values over time. The latter is explained in more detail below.

Figure 5 Quantile–quantile (QQ) plots for different numbers of input observations. For example, 3 shows the quantiles of all SVODI values where three input indices contributed. Any shows the quantiles for any number of observations. N(0,1), the diagonal, shows the theoretical normal distribution for comparison. (a) The QQ plot of SVODI. (b) For comparison the quantiles of an index generated by simple averaging.

3.2.1 Temporal shift determination

Is it of interest whether events can be seen first in the microwave or optical domain and how large their temporal difference is for a variety of applications, such as drought prediction. For this purpose the temporal shifts between SVODI and VHI and soil moisture anomalies are determined by finding the temporal lag at which a dataset pair correlates most strongly. This is done by grid search, calculating the correlation coefficient for every shift within a ± 8-week window and selecting for each location the shift with the highest correlation. An 8-week window was chosen because the shifts were found to be almost always within that, even when searching in a larger window. Then, all results are filtered for unreliable results: all locations where the correlation as a function of temporal shift exhibits multiple local maxima, detected by counting the number of sign changes of the first derivative, or where the maximum correlation coefficient is less than zero are masked out.

Figure 6 Spatial correlation coefficient of SVODI vs. VHI over time (a) and per season (b), based on weekly data.

3.2.2 Extreme values over time

The plots showing the frequency of extreme values over time are inspired by the plots in Van Der Schrier et al. (2013). With these plots it is possible to visualize whether the temporal abundance of extreme values is similar for the different datasets.

For visualization, the data are first standardized to N(0,1) using Eq. (1). For example, VHI, TCI and VCI are originally scaled from 0 to 100, and even the N(0,1) distributed indices (e.g. SVODI) are sometimes temporally downsampled, which leads to a reduced variance and therefore requires re-standardization. The percentage of extreme values is then calculated as follows: for a geographical region, for each time step, all pixels with a value greater than 1 and greater than 2 are counted and divided by the total available data points available for that time step.

4 Results and discussion

4.1 Technical analysis

4.1.1 Illustration of the methods by means of an example time series

Figure 3 shows the creation of SVODI at various steps for an exemplary location. The original series have a visible bias between values of the same band (panel a), which is corrected with the CDF matching (panel b). Note that the scaling is done individually for each frequency band and not to a common reference used for all together. For example, AMSR2, WindSat and TMI Ku-band observations are matched to AMSR-E Ku-band observations. Then, an index is created for each sensor and band, and the multiple indices are finally merged into SVODI (panel c). The prior distribution of SVODI is unaffected by the number of observations available at a certain date as the number of available observations varies from two to eight (panel d), but the SVODI time series shows no breaks between periods with different sensor availability.

4.1.2 Extreme values over time

Prior to the SVODI calculation, all datasets are detrended. On a global scale it is therefore expected that the percentage of extreme SVODI values, both positive and negative, is more or less constant over time. Indeed, there seems to be no drastic systematic increase or decrease in the percentage of pixels with |SVODI| greater than 1 or 2 (Fig. 4a). This indicates that even though the number of sensors contributing to SVODI changes over time, this does not lead to massively more or fewer extreme events in a given time period. In contrast, Fig. 4b shows the percentage of extreme values if SVODI were generated by simply averaging the individual indices per band and sensor and standardizing the average again. In this case, the periods with few sensors (pre-2002, after 2017) are much more likely to be extreme than the periods with more sensors. This shows that our probabilistic merging method is necessary to compare the frequency of extreme events across different periods.

4.1.3 Impact of the number of input datasets on SVODI

If, compared to a simple VODCA standardization, the number of input sensors were to have no effect on the SVODI distribution, then this should be standard normally distributed, irrespective of the number of input sensors. Figure 5a shows the quantiles of SVODI with respect to the quantiles of a standard normal distribution for different numbers of input sensors. Note that this figure is based on a random 20 % of all data as the whole dataset is too large to be loaded at once. Also, each SVODI value can only be part of one group. Hence, each group distribution is computed from values from different dates.

Generally, SVODI is normally distributed regardless of the number of input sensors used for its computation. Only for extremely low values is a small difference observed. Very low values that are also the result of many sensors have a slight positive bias, while for very high values this discrepancy does not occur. The cause for this problem is not fully understood, but it is assumed to be related to numerical instability of very low values as p for days with many observations can become very low. As a simplified numerical example, if a SVODI value is the result of eight individual uncorrelated indices and all of them are 0 (indicating average vegetation conditions for all sensors), then the resulting p is CDF$(\mathrm{0}{)}^{\mathrm{8}}={\mathrm{0.5}}^{\mathrm{8}}=\mathrm{0.004}$. This very low value has then to be scaled back to 0.5, and therefore even minor deviations in p can lead to a substantially different p_scaled. This example shows that the computation can involve very low values and can therefore become unstable.

Figure 5b shows the quantiles if, rather than our probabilistic merging method, the simple mean of each individual index were used and the result were standardized. It shows that the values of this aggregate index would strongly depend on the number of input observations as it is much more likely to be extreme if only a few input observations are available. This shows that SVODI's dependency on the number of observations is almost removed and a lot lower than if one were to use a simple standardization.

Figure 7 Correlation coefficient without temporal shift between SVODI and VCI (a), TCI (c) and VHI (e) and the temporal shift at which a maximum correlation is obtained (b, d, f). All results are based on weekly means; positive (green) values in the shift plots indicate that anomalies are visible in VCI, TCI or VHI before SVODI.

Figure 8 Correlation (a, c, e, g) and temporal shift (b, d, f, h) between SVODI and soil moisture anomalies from ERA5, based on weekly data. Positive (green) values in the shift plots indicate that anomalies are earlier visible in soil moisture than in SVODI.

4.1.4 Quality change over time

Of interest is whether the quality of SVODI changes over time. There exist no ground-based validation data for VOD; therefore a direct validation with some reference data is not possible. Instead, the spatial correlation to VHI over time is used as an indicator of whether SVODI is performing differently during different periods. For each time step, the correlation between the global SVODI and VHI images is calculated. If the signal-to-noise ratio of any of the two datasets were to change, so would the correlation between them.

The spatial correlation varies quite strongly over time (Fig. 6a). This variance is likely due to the varying number of events over time. There is a statistically significant positive trend. This is expected as more and better sensors are progressively becoming available in later years for both SVODI and VHI, resulting in higher-quality data for both. In the future, as more advanced sensors are launched into orbit, the two datasets are expected to further converge.

There is no apparent seasonal dependency of the correlation (Fig. 6b), indicating that the increase in correlation is not an artefact of long-term seasonal surface changes.

Figure 9 Mean correlation coefficient if no temporal shifting is done (a), maximum correlation obtained by temporal shifting (b) and the corresponding temporal shift in weeks between SVODI and soil moisture anomalies (c) at different ERA5 depths. The depths are swvl1 0–7, swvl2 7–28, swvl3 28–100 and swvl4 100–289 cm. Panel (a) also shows the correlation with the surface soil temperature, stl1.

Figure 10 Percentage area of SVODI greater/smaller than 1/−1 and 2/−2 for central Australia together with the SOI (a), SVODI (b, d) and standardized precipitation anomalies (c, e) for November 2010 and December 2019.

4.2 Data analysis

4.2.1 Comparison of SVODI and VHI

SVODI is compared to VCI; TCI; and their composite, VHI, to explore their similarities and differences. By comparing SVODI to all three, it is possible to evaluate how it relates to the impact of stress on either vegetation “greenness” (VCI) or temperature (TCI) or a combination of both (VHI).

VCI and SVODI correlate quite strongly, especially in semi-arid climates (Fig. 7a). This is in line with previous studies comparing VOD anomalies with LAI, which is also derived from optical data (Moesinger et al., 2020; Jones et al., 2011). The pattern is at least partially due to the more distinct inter- and intra-annual variability in vegetation activity in semi-arid regions. Vegetation in semi-arid regions is highly susceptible to increased or decreased precipitation, and as such these regions experience stronger anomalies than regions where water is abundant, leading to a higher signal-to-noise ratio. This then leads to higher correlations between the two indices.

TCI and SVODI correlate positively in semi-arid regions (Fig. 7c) because high temperatures lead to unfavourable vegetation conditions and vice versa for low temperatures. Likewise they correlate negatively in cold regions where positive temperature anomalies lead to favourable vegetation conditions and vice versa for low-temperature anomalies. (Papagiannopoulou et al., 2017). SVODI correlates more strongly with VCI than with TCI. This makes sense as both SVODI and VCI represent changes in biogeophysical properties of the canopy as a result of anomalous environmental conditions, whereas TCI mirrors a cause of these changes.

Figure 11 Percentage area of SVODI (a), standardized VHI (b) and scPDSI (c) greater/smaller than 1/−1 and 2/−2, respectively, over time for the northern South America AR6 region (Iturbide et al., 2020). All datasets are downsampled to the monthly resolution of scPDSI.

Since VHI is the average of TCI and VCI, SVODI also correlates more strongly with VCI than with VHI (Fig. 7e), especially in cold regions.

Anomalies generally occur first in TCI followed by VCI and SVODI (Fig. 7b, d, f), which is expected as there is a causal relationship between prolonged high temperatures in semi-arid regions and a subsequent reduction in vegetation health and vice versa for low temperatures. VCI generally leading SVODI indicates that changes in greenness occur before changes in vegetation water content. This is to be expected as reduced photosynthesis is one of the first responses of plants to heat stress, in part due to decay of photosynthetic pigments (Larcher, 2000; Zhao et al., 2020). The plant water content drops more slowly as plants are able to regulate the rate of transpiration and respiration to balance water loss under transient or mild heat stress (Zhao et al., 2020).

4.2.2 Relation between SVODI and soil moisture at different depths

Correlations between SVODI and soil moisture anomalies at various depths were calculated to determine the connection between the different depths and the vegetation condition. SVODI and upper level soil moisture anomalies (0–7 and 7–28 cm) correlate most strongly with each other in areas where vegetation growth is limited by water availability (Fig. 8a, c). This is in line with a previous study comparing soil moisture to VOD anomalies (Konings et al., 2021). In areas where vegetation growth is limited by temperature (e.g. see Hashimoto et al., 2019; Papagiannopoulou et al., 2017), there is generally a negative correlation between SVODI and soil moisture at all levels. This makes sense, as increased precipitation, increased soil moisture and increased cloud coverage and therefore lower temperatures are all linked to each other in these regions. Additionally, increased soil moisture leads to decreased LST by controlling the partitioning between sensible and latent heat fluxes (Huang and van Den Dool, 1993). The lower LST then leads to less favourable vegetation conditions in cold regions. This can also be seen if the mean correlation per land cover and depth is considered (Fig. 9a). In land cover types that are generally found in cold regions, such as needleleaf forests, surface soil moisture anomalies correlate negatively with SVODI and slightly positively with temperature anomalies. The opposite is the case for land cover types generally found in warmer climates (broadleaf forests, shrubs), where more water and lower temperatures lead to more favourable vegetation conditions. Roughly the same pattern is also visible if the analysis is repeated for each season separately (DJF, MAM, JJA, SON¹), but the pattern shifts with the seasons (results not shown). For example, in Europe, SVODI correlates negatively with soil moisture in DJF due to the temperature being the limiting factor but positively in JJA when water is the limiting factor.

Correlation coefficients between SVODI and soil moisture anomalies decrease with increasing depth when no lag optimization is performed (Fig. 9a). However, by optimizing the lag, the highest correlation coefficients are obtained for 7–28 and 28–100 cm soil moisture (Fig. 9b), even though overall the differences are very small. This agrees with a study showing 7–28 cm to be the most relevant water reservoir for vegetation productivity, particularly in semi-arid regions (Li et al., 2021). The optimized correlation coefficients are highest for land cover types typically found in arid regions, such as sparse vegetation and shrubs, while they are lowest in needleleaf forests, which are more often found in regions where temperature is usually more limiting than water.

There is a clear relationship between soil depth and temporal shift (Fig. 9c), with soil moisture anomalies in layer 3 generally preceding SVODI and soil moisture anomalies in layer 4 generally following SVODI. This is mostly likely due to bottom-level soil moisture levels lagging behind the top ones, which makes sense if moisture is modelled as a “bucket” model where the top layers are filled and depleted first before the lower layers.

4.2.3 Sensitivity of SVODI to Australian interannual precipitation variability

The Australian summers of 2010 and 2019 were marked by exceptionally high and low precipitation (Fig. 10c and e), respectively. The year 2019 was also exceptionally hot, and widespread wildfires occurred (Dunn et al., 2020). As the Australian vegetation is strongly limited by water availability, one would expect the vegetation moisture content in these years to react similarly. Indeed, the SVODI of 2010 and 2019 was exceptionally high and low, respectively (Fig. 10a), and the spatial patterns of SVODI and precipitation anomalies are very similar in these periods (Fig. 10b and d). This illustrates that SVODI in this water-limited region behaves as expected and is useful to monitor the state of the vegetation.

4.2.4 Effects of drought on vegetation in the Amazon

In the Amazonian rainforest, the extreme values of SVODI, VHI and scPDSI do not agree with each other (Fig. 11). SVODI seems to have more high-frequency variability but less low-frequency variability than VHI and scPDSI. One possible explanation might be the originally higher spatial resolution of the VHI and scPDSI products, which, after downsampling to the spatial resolution of SVODI, lead to very low noise levels in these products. Another possibility is that SVODI reacts more quickly to surface changes than the other two indices.

The effects of droughts in the Amazon forest is a highly discussed topic and very challenging (Samanta et al., 2012). Highly discussed droughts occurred in 2005, 2010 and 2015 (Panisset et al., 2018; Liu et al., 2018; Lewis et al., 2011; Samanta et al., 2010; Janssen et al., 2021). Some studies argue that the tropical vegetation is light-limited, and as such a decrease in precipitation, which corresponds to a decrease in cloud cover, leads to greening during droughts (Huete et al., 2006; Saleska et al., 2007). Others argue that the greening is due to artefacts from atmospheric effects and changing sun-sensor geometry (Samanta et al., 2010; Morton et al., 2014). Our results do not give a definitive answer to this discussion as the patterns are very ambiguous.

4.2.5 Climate oscillation relations to SVODI

The correlation between SVODI, VHI and scPDSI and the Southern Oscillation index (SOI) (Fig. 12, left column) and dipole mode index (DMI) (Fig. 12, right column) is calculated to study how global vegetation is affected by tropical climate oscillations.

Figure 12 Correlation coefficients between SVODI, VHI and scPDSI (top to bottom) vs. SOI (a, c, e) and DMI (b, d, f).

SVODI and VHI correlations show a similar pattern. The highest correlation coefficients to SOI are found in eastern Australia, where vegetation is heavily influenced by the El Niño–Southern Oscillation (ENSO) (Liu et al., 2009; Miralles et al., 2014), with high SOI values being linked to increased precipitation. While scPDSI also agrees with VHI and scPDSI in Australia, it also has a high positive correlation in northern South America and a negative correlation in southern North America and western Asia.

The correlations between SVODI and VHI and DMI are similar and show the greatest magnitude in southern Australia, while in most other regions they are close to zero. Correlations between scPDSI and DMI show a more distinguished spatial pattern, including strong positive correlations in north-western Russia. In this region, vegetation growth is limited by temperature and not precipitation, and for this reason no corresponding patterns can be found in the SVODI or VHI plots.

5 Conclusions

SVODI is a microwave-based vegetation condition index that shows similar patterns to existing optical indices and follows soil moisture in semi-arid regions. It extends the current range of available remote sensing datasets that allow the observation of anomalous vegetation states and increases our understanding of global vegetation dynamics. SVODI patterns are reasonable compared to patterns of VHI, TCI, VCI and soil moisture, but anomalies occur later, which might be an issue for near-real-time applications.

With the exception of extreme low values, the proposed index generation method works well at combining different indices. The merging method itself is not limited to VOD and can potentially be applied to combine arbitrary normally distributed indices. Therefore, this method might find applications in various disciplines. Further efforts will focus on increased numerical stability of the calculations and updating SVODI with more recent observations.