Contribution of Land Surface Temperature ( TCI ) to Vegetation Health Index : A Comparative Study Using Clear Sky and All-Weather Climate Data Records

The Vegetation Health Index (VHI) is widely used for monitoring drought using satellite data. VHI depends on vegetation state and thermal stress, respectively assessed via (i) the Vegetation Condition Index (VCI) that usually relies on information from the visible and near infra-red parts of the spectrum (in the form of Normalized Difference Vegetation Index, NDVI); and (ii) the Thermal Condition Index (TCI), based on top of atmosphere thermal infrared (TIR) brightness temperature or on TIR-derived Land Surface Temperature (LST). VHI is then estimated as a weighted average of VCI and TCI. However, the optimum weights of the two components are usually not known and VHI is usually estimated attributing a weight of 0.5 to both. Using a previously developed methodology for the Euro-Mediterranean region, we show that the multi-scalar drought index (SPEI) may be used to obtain optimal weights for VCI and TCI over the area covered by Meteosat satellites that includes Africa, Europe, and part of South America. The procedure is applied using clear-sky Meteosat Climate Data Records (CDRs) and all-sky LST derived by combining satellite and reanalysis data. Results obtained present a coherent spatial distribution of VCI and TCI weights when estimated using clearand all-sky LST. This study paves the way for the development of a future VHI near-real time operational product for drought monitoring based on information from Meteosat satellites.


Introduction
Climate Data Records (CDRs) [1] have been continuously and incrementally used in the last decade, and agencies such as the European Space Agency (ESA) or the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT) provide datasets [2][3][4] for a diverse range of Essential Climate Variables (ECVs) that can be generated using satellite observations.These ECV CDRs cover parameters in various domains of the climate system, like the atmosphere (e.g., total column water vapor [5], Top-of-Atmosphere (TOA) radiance [6]), the ocean [7], and the land (e.g., albedo [8,9], Land Surface Temperature (LST) [10]).LST has recently been added as an ECV following the recognition of its role in processes within the interface between the soil and the atmosphere.The use of ECV CDRs-consisting of land surface and/or atmospheric variables over large regions with spatially and temporally homogeneous information over decades of Earth observation-may contribute to a more performant monitoring of climate elements like snow cover and extreme events like drought episodes [11][12][13].Indeed, the latter are a recurrent climate event with severe impacts and consequences on society [14,15] that may result in famine, death [16], and large economic losses [17].Several different drought indices were developed in the last decades with the aim of characterizing and monitoring the different types of drought [18][19][20].Examples of such indices are the Palmer Drought Severity Index (PDSI) [21], the Standardized Precipitation Index (SPI) [22], and the Standardized Precipitation-Evapotranspiration Index (SPEI) [23,24].Satellite-based remote sensing is also largely suited for monitoring drought conditions since its instruments and sensors allow for both regional and global scale observation of the land surface with spatially continuous measurements and a high repeat visit frequency in the case of geostationary satellites.Several indices based on satellite thermal infrared (TIR) data were developed through the years, such as the Apparent Thermal Inertia (ATI) [25,26], the Temperature Vegetation Dryness Index (TVDI) [27], or the Vegetation Health Index (VHI) [28,29], among others [30].A more performant drought monitoring at the continental or global scales will contribute to the adoption of more efficient measures to mitigate the adverse impacts.In this context, a crucial role is played by indices that, although suitable to be applied at the global scale, are tuned for the different biomes and regions so that they optimally reflect the stress of vegetation.Such tuning will benefit from CDRs since they provide consistent information at the large scale, whereas traditional meteorological observations are usually sparser and more limited.
VHI is a widely used remote sensing-based drought index designed as the weighted sum of two components: the Vegetation Condition Index (VCI) and the Thermal Condition Index (TCI).The first component characterizes moisture conditions and is typically based on information from the visible and near infra-red windows of the electromagnetic spectrum, whereas the latter characterizes the thermal condition and is based on information from the thermal infra-red window.The Normalized Difference Vegetation Index (NDVI) and LST (or TOA brightness temperature) are commonly used to estimate VCI and TCI, respectively.The optimum weights that are to be attributed to VCI and TCI are usually not known and VHI is typically estimated by assuming equal weights of 0.5 to both components.Results from a recent study [13] point out that it is possible to improve the estimation of VCI and TCI weights by comparing the resulting VHI against a non-satellite derived multi-scalar drought indicator such as SPEI.The study was based on the use of an LST CDR developed by Princeton University (PU hereafter) [31,32] where satellite observations were merged with reanalyzes to derive an all-weather LST, assuming land surface emissivity equal to one.However, most satellite LST products are based on infra-red measurements and are therefore available under clear-sky conditions only.LST datasets are also likely to vary in terms of algorithms and input data or assumptions.The aim of this work is to analyze the sensitivity of drought assessment, via an improved VHI index as the one derived in Ref. [13], to the type LST dataset considered, namely the all-weather or the clear-sky ones.This is particularly relevant, as the use of a methodology based on VHI, and therefore on LST and on a vegetation index such as NDVI, is particularly suitable for the monitoring of drought over large areas using remote sensing data alone.For that purpose, we compare VHI derived with different LST datasets, namely the PU LST and the two clear-sky LST CDRs recently released by the EUMETSAT Satellite Application Facility on Climate Monitoring (CM-SAF), both derived from Meteosat First (MFG) and Second (MSG) Generations [10,33].
The Meteosat disk encompasses the African continent, southern Europe, and the northeast of Brazil.The Sahara Desert was not included in the analysis of drought conditions because of its persistent very arid character; in turn, the tropical rain forests and moister woody savannas found in equatorial central Africa are usually affected by very cloudy conditions and therefore fewer and less accurate satellite observations of the surface are available.Most of the remaining regions are however plagued by recurrent drought episodes: Iberia [34,35], Euro-and African-Mediterranean region [36][37][38], eastern, southern, and the Sahel regions of Africa [39], and north-eastern Brazil [40,41].These are, consequently, regions of interest for this kind of analysis.Furthermore, drought frequency and severity are expected to increase due to climate change in regions like southern Europe and the Sahel [42-47], making the region encompassed by the Meteosat disk a hotspot for future drought monitoring.
The main objectives of this study are twofold: (1) to use the previously developed methodology to estimate contributions of thermal and vegetation stress to vegetation health [13] and (2) to understand whether, by using different LST CDRs-while preserving NDVI and SPEI-there are large differences from the optimum design of VHI, namely, if different weights should be attributed to the contributions of thermal stress and vegetation state in order to obtain an optimal vegetation health index.The behavior of contributions of thermal and vegetation stress for different biomes within the Meteosat disk is analyzed and three severe drought episodes, two in Iberia and one in Morocco are studied.This analysis may be viewed as a precursor for the development of a real-time product for drought monitoring using MSG and the future Meteosat Third Generation (MTG) instrument [48,49].The paper is organized as follows: in Section 2, the CDRs used in this study are presented and the methodology to estimate the VHI contributions is described; in Section 3, the optimum VCI/TCI weights for the computation of VHI derived using three different LST CDRs are intercompared and a biome analysis to assess dependencies of the contributions is performed; Section 4 includes the discussion of the results; and finally, Section 5 presents the overall conclusions and possible future work.

Study Area and Biome Classification
The geographical region of interest for this study covers the Meteosat viewing disk from −65 • S to 65 • N in latitude and −65 • W to 65 • E in longitude (Figure 1).
The study is performed for land surfaces within the African continent, southern Europe, part of the middle East, and north-eastern Brazil.As mentioned above, one of the objectives of this study is to assess variations of the contributions of thermal stress and vegetation state to VHI over different biomes.The land biome classification used here is a variant of the 2006 GlobCover product [50] developed for the Advanced Along-Track Scanning Radiometer (AATSR), with a spatial resolution of 1/120 • re-gridded from the original 1/360 • GlobCover resolution [51].This land cover dataset is stratified into a total of 27 land and inland water/coastal water biomes and one open ocean biome and takes into consideration changes in emissivity among biomes.

LST, NDVI Climate Data Records, and SPEI Dataset
Here, three LST CDRs are considered.The first LST dataset was developed at Princeton University (PU) [31,32] and consists of hourly data defined on a global high-resolution grid (0.5 • ) for a period spanning from January 1979 to December 2009.This all-sky LST dataset relies on information from the High-Resolution Infrared Radiation Sounder (HIRS) retrieved LST merged through a Bayesian postprocessing methodology with the National Centers for Environmental Prediction (NCEP) Climate Forecast System Reanalysis (CFSR) LST with the aim of filling unobserved estimates of satellite retrievals.The dataset was comprehensively validated with in situ LST observations [31].In order to provide appropriate consistency between datasets, the authors [31] used a unitary emissivity to retrieve satellite and reanalyzed LSTs.The merged variable therefore represents a surface brightness temperature.However, the term LST that is used for this variable by the dataset authors is kept here.
The remaining two LST CDRs used in this work are developed at CM-SAF (in a joint effort with the EUMETSAT SAF on Land Surface Analysis, LSA-SAF) and are derived from the geostationary MFG and MSG satellites using the so-called Physical Mono-Window (PMW) and Statistical Mono-Window (SMW) algorithms [10,33].The first algorithm is based on radiative transfer simulations inverting the radiative transfer equation and Planck's distribution [10,52,53], whereas the latter relies on the linear relationship between LST and top of atmosphere brightness temperatures and consists on linearizing the radiative transfer equation while maintaining an explicit dependency on surface emissivity [10].Both datasets rely on observations of a single channel in the thermal infra-red window, so that the same algorithm can be used for both MFG and MSG.The two CDRs are provided in a regular latitude and longitude grid with spatial resolution of 0.05 • with the geographical coverage detailed in Section 2.1.CM-SAF PMW and SMW LST CDRs cover the period ranging from January 1991 to December 2015, and the monthly diurnal cycle is used in this work.These datasets consist of satellite data alone, and therefore it is worth stressing that they are clear-sky estimates and that, in this case, a map of variable emissivities [54] is used.The CM-SAF LST products were compared [55] in terms of: (1) precision and accuracy against more than 50,000 in-situ LST measurements and showed errors below 1 K for the monthly CDRs; and (2) decadal stability against ERA-Interim skin temperatures and the Moderate Resolution Imaging Spectroradiometer (MODIS) LST where errors were below 0.8 K.
The Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI used in this study was processed and produced by the NASA/Goddard Space Flight Center's Global Inventory Modelling and Mapping Studies (GIMMS) Group [56].This NDVI CDR, consisting of 8 km spatial resolution and bi-monthly temporal resolution values, starts in January 1982 and is the longest and most complete NDVI record available.
The multi-scalar drought index SPEI used in this study, from the Consejo Superior de Investigaciones Científicas (CSIC, Spain; http://spei.csic.es/index.html)[23,24,57], relies on Climate Research Unit (CRU) TS 3.23 datasets of temperature and precipitation [58] and is based on a normalization of the simple water balance developed by Thornthwaite [59].This monthly global dataset covers the period ranging from January 1901 to December 2015 with a spatial resolution of 0.5 • , and as a multi-scale character with available time-scales between 1 and 48 months.
In order for the datasets to be consistent, the temporal coverage is truncated to span the period between January 1991 and December 2009 for LST and NDVI datasets, and the period between December 1990 and November 2009 for the SPEI dataset (see Section 2.2 for details about this 1-month lag).Furthermore, PU LST is aggregated into monthly diurnal cycle composites to match the CM-SAF LST datasets.Monthly means of the three LST datasets are then estimated by averaging the month diurnal cycle.NDVI bi-monthly values are averaged to get monthly values.NDVI and CM-SAF LSTs are interpolated to a 0.5 • × 0.5 • spatial grid.Finally, SPEI time-scales of 1-, 3-, 6-, 9-, and 12-months are chosen.Linear trends on NDVI associated to the region's greening or browning [60,61], and those associated to global warming on LST, are removed by applying a linear regression algorithm.Furthermore, the biome classification is also re-gridded to a 0.5 • × 0.5 • spatial resolution by choosing the dominant biome.
In summary, this study considers four monthly Climate Data Records (three LST and one NDVI) and one SPEI dataset collocated over a grid of 0.5 • and spanning a period of 19 years.

Vegetation Health Index
VHI characterizes the health of the vegetation by assuming that stressed conditions are linked to lower than normal NDVI and higher than normal temperature [28,29]; it may be written as: where VCI is the Vegetation Condition Index, TCI is the Thermal Condition Index and α is a parameter that quantifies the contribution of each component to the overall vegetation health.VCI and TCI are defined as: where NDVI, NDVI min , and NDVI max are the value for a specific month and pixel and the minimum/maximum value for the same month and pixel over the climatological period of study, respectively.The same notation is valid for LST, LST min , and LST max .
Here, three datasets of TCI are estimated: with PU LST, CM-SAF PMW LST, and CM-SAF SMW LST; while VCI is estimated with the GIMMS NDVI.Consequently, three datasets of VHI are also estimated.
Parameter α in Equation ( 1) takes a value between 0 and 1.Since normally there is no a priori knowledge of the actual contribution of the moisture and temperature conditions to the health of the vegetation in a given region, the value of α is generally taken as 0.5.However, as pointed out in Ref. [13], with the help of a multi-scalar drought estimator, such as SPEI, it is possible to estimate an optimal value of α for the different biomes and therefore assess the relative contributions of VCI and TCI for different regions.In the next section, the methodology to estimate α is presented.

Estimating Moisture (VCI) and Temperature (TCI) Contributions to VHI
The procedure to estimate parameter α in Equation (1) follows the one described in Ref. [13].
Here, VHI is first estimated with the three different LST datasets by assuming the traditional α = 0.5, i.e., equal contributions of the moisture (VCI) and temperature (TCI) conditions to vegetation health.For each pixel, the month of climatological maximum NDVI is chosen, and time series of VHI of that month are correlated with time series of SPEI from the previous month at time-scales of 1-, 3-, 6-, 9-, and 12-months.The month of maximum NDVI is chosen since it is a good indicator to assess the influence of droughts measured with SPEI in the growing season to the vegetation health of the peak vegetation productivity period.The choice of this month also simplifies the analysis in terms of the amount of information to process.Then, the SPEI time-scale that presents the maximum correlation with VHI is chosen [62] and a p-value of 0.1 for the two-tailed t-test is used to reject the hypothesis of no correlation.
For pixels where the obtained correlation is significant, an iterative maximization of the correlation between VHI and SPEI is performed, using values of α ranging from 0.025 to 0.975 in steps of 0.05.This may be regarded as an adjustment of the contributions of VCI and TCI for pixels where the traditional definition of VHI is positively correlated with SPEI. Figure 1 shows the resulting maps of α and respective SPEI time-scale.correlation with VHI is chosen [62] and a p-value of 0.1 for the two-tailed t-test is used to reject the hypothesis of no correlation.
For pixels where the obtained correlation is significant, an iterative maximization of the correlation between VHI and SPEI is performed, using values of α ranging from 0.025 to 0.975 in steps of 0.05.This may be regarded as an adjustment of the contributions of VCI and TCI for pixels where the traditional definition of VHI is positively correlated with SPEI. Figure 1 shows the resulting maps of α and respective SPEI time-scale.

Comparison of α Estimated with the Different LST CDRs
The optimal values of parameter α estimated using the three LST datasets are shown in Figure 1.These correspond to the estimate of VHI that leads to the best correlation with the relevant SPEI (scale also shown in Figure 1).Only pixels with statistically significant correlations (p < 0.1) are shown.These are confined to Mediterranean Europe, Southern Africa, North-East Brazil, and a few scattered regions in the Sahel and sub-Saharan Africa.As already pointed out, areas characterized by desert/very arid climates (Sahara, Arabian, and Namib Deserts and the eastern part of the horn of Africa) were previously masked out.Forests in central Africa and parts of Europe are characterized by non-significant correlations between VHI and SPEI (Figure 1) [13].
The patterns of both α and SPEI time-scale are very similar among the three datasets.Most of the pixels show high values of α, while lower values are confined to relatively small clustered regions (e.g., in east and south Africa and East and Northeast Brazil).Concerning SPEI time-scales, the lower range of 1-and 3-month dominates over most of the Meteosat disk.Nevertheless, longer time-scales are organized in small clusters scattered throughout the region.It is worth noting that these features are present in the maps derived from the three LST datasets.
The magnitude and spatial distribution of α are coherent among these datasets, when either the LST datasets were derived from the same observations with different algorithms (the case of CM-SAF LST CDRs) or from different observations and methodologies (PU and CM-SAF LST CDRs).
For a closer look, Figure 2 presents a comparison among values of parameter α as estimated with the three LST CDRs (PU LST, CM-SAF PMW LST, and CM-SAF SMW LST).The left panels compare the distributions of α estimated with CM-SAF PMW LST versus PU LST (top), with CM-SAF SMW LST versus PU LST (middle), and CM-SAF PMW versus CM-SAF SMW (bottom).Right panels show histograms of the differences between the same pairs of variables.Results point out that values of α are well distributed along the diagonal (solid black line) in all three cases.Furthermore, the agreement between α estimated with CM-SAF datasets is higher, as expected, since these correspond to LST derived from the same observations-both clear-sky only.
Values of Root Mean Squared Deviation (RMSD) of α derived with PU LST and clear-sky PMW and SMW LSTs are of about 0.2.Furthermore, it is worth noticing that the differences between α estimated with CM-SAF LSTs are smaller, with a RMSD of about 0.1.The average differences of α estimated with the different LST datasets are negligible, meaning that the use of these datasets is not introducing any overall systematic error in this type of analysis.
Figure 3 shows the differences ∆ρ between correlations of VHI and SPEI when α is estimated with the iterative maximisation (ρ it ) method and when α is taken as 0.5 (ρ 0.5 ): Since the iterative procedure is designed to obtain the maximum value of correlation, ∆ρ is always greater to or equal than 0. It is worth noting that regions where ∆ρ is larger match those where values of α depart more from 0.5.Correlation differences present very similar standard deviation (about 0.06) and bias (about 0.05) for all datasets.Nevertheless, when using PU LST, the patterns over South Africa present larger increases in the correlations, and therefore a better VHI optimization, than CM-SAF LSTs (the histogram also shows less pixels in the first bin).

Biome Stratification
Results presented in the previous section indicate that α is organized in clusters of similar values over the different regions.With the aim of assessing whether these clusters are associated with land cover, the ATSR LST biome classification derived from GlobCover was used to analyze the distribution of land cover among pixels with significant correlation between VHI and SPEI (Figure 4).From the 27 biomes of the original classification, 21 are found among the significant pixels and only 9 contain more than 100 pixels (Figure 4).These 9 biomes are therefore selected to assess the dependence of VHI on the moisture and temperature conditions and to assess the differences between the contributions of these two components when different LST CDRs are used to estimate TCI.
Figure 5 presents results obtained for the biome stratification in the form of box plots for α (top), Δρ (middle) and ΔVHI (bottom) as derived using PU LST (black), CM-SAF PMW LST (blue) and CM-SAF SMW LST (red) to estimate TCI.As for Δρ, ΔVHI represents the differences between VHI using α estimated with the iterative maximization based on SPEI (VHI ) and VHI using α set to the standard value of 0.5 as originally proposed by Kogan (VHI ): The box delimits the 25th and 75th percentiles and the inner segment represents the median, whereas the whiskers represent percentiles 10 and 90.Table 1 presents percentiles 25 and 75 and the median of the distribution of α estimated with the three LST datasets, for the different biomes presented in Figure 4.

Biome Stratification
Results presented in the previous section indicate that α is organized in clusters of similar values over the different regions.With the aim of assessing whether these clusters are associated with land cover, the ATSR LST biome classification derived from GlobCover was used to analyze the distribution of land cover among pixels with significant correlation between VHI and SPEI (Figure 4).From the 27 biomes of the original classification, 21 are found among the significant pixels and only 9 contain more than 100 pixels (Figure 4).These 9 biomes are therefore selected to assess the dependence of VHI on the moisture and temperature conditions and to assess the differences between the contributions of these two components when different LST CDRs are used to estimate TCI.
Figure 5 presents results obtained for the biome stratification in the form of box plots for α (top), ∆ρ (middle) and ∆ VHI (bottom) as derived using PU LST (black), CM-SAF PMW LST (blue) and CM-SAF SMW LST (red) to estimate TCI.As for ∆ρ, ∆ VHI represents the differences between VHI using α estimated with the iterative maximization based on SPEI (VHI SPEI−based ) and VHI using α set to the standard value of 0.5 as originally proposed by Kogan (VHI K ): The box delimits the 25th and 75th percentiles and the inner segment represents the median, whereas the whiskers represent percentiles 10 and 90.Table 1 presents percentiles 25 and 75 and the median of the distribution of α estimated with the three LST datasets, for the different biomes presented in Figure 4.As shown in Table 1 and as suggested by the box plots of Figure 5 (top panel), the distributions tend to be shifted towards larger values of α (box plots on the top panel of Figure 5; Table 1) with the median being above 0.5 for most cases.Over open broadleaved deciduous forests (biome 07) with α estimated with PU LST and over closed to open shrublands (biome 13) with α estimated with both CM-SAF LSTs, the median of optimum α falls slightly below 0.5 (and is slightly above 0.5 for the same biomes but when computed with a different LST dataset).The highest values of α may be found in southwestern Africa over closed to open grasslands (biome 14) as well as over Iberia and northern Africa over sparse vegetation (biome 15, which also contains grasslands, as well as shrubs and woody vegetation).The closed to open grasslands present medians of α near 0.9 and percentile 25 is higher than 0.6 for the three LST datasets and, when α is estimated with PU LST percentile 10 is also above 0.5, whereas when α is estimated with CM-SAF LSTs it is close to 0.4.Sparse vegetation shows medians of α of about 0.7 for the three LST datasets (with PU LST showing α higher by 0.05 than the remaining).As shown in Table 1 and as suggested by the box plots of Figure 5 (top panel), the distributions tend to be shifted towards larger values of α (box plots on the top panel of Figure 5; Table 1) with the median being above 0.5 for most cases.Over open broadleaved deciduous forests (biome 07) with α estimated with PU LST and over closed to open shrublands (biome 13) with α estimated with both CM-SAF LSTs, the median of optimum α falls slightly below 0.5 (and is slightly above 0.5 for the same biomes but when computed with a different LST dataset).The highest values of α may be found in southwestern Africa over closed to open grasslands (biome 14) as well as over Iberia and northern Africa over sparse vegetation (biome 15, which also contains grasslands, as well as shrubs and woody vegetation).The closed to open grasslands present medians of α near 0.9 and percentile 25 is higher than 0.6 for the three LST datasets and, when α is estimated with PU LST percentile 10 is also above 0.5, whereas when α is estimated with CM-SAF LSTs it is close to 0.4.Sparse vegetation shows medians of α of about 0.7 for the three LST datasets (with PU LST showing α higher by 0.05 than the remaining).Table 1.First quartile (p25), median (p50), and third quartile (p75) of the distribution of α as estimated with PU LST (PU), CM-SAF PMW LST (PMW), and CM-SAF SMW LST (SMW) for the different biomes (identified by the labels described in Figure 4).Differences in the distribution of α are also found when comparing results for the three LST datasets (Figure 5, top; Table 1): ( 1 On the other hand, the distribution of ∆ρ (Figure 5, bottom panel) shows larger correlation increases for biomes characterized by highest values of α, namely grasslands and sparse vegetation.For the latter biomes, results with PU LST show larger increases in the correlation than those obtained with CM-SAF LSTs.Consequently, the greater differences in VHI are also found over these biomes and for PU LST.

Comparison of VHI for Different Case Studies
Three cases of documented [63][64][65] severe drought episodes over Iberia and Morocco are studied: (1) the winter drought over southern Iberia in 1998/1999; (2) the spring drought over Iberia in 2005; and (3) the spring drought of 2007 in Morocco.For this purpose, two areas inside Iberia and Morocco were chosen.Figure 6 shows the time-series of VHI estimated with the traditional rationale using equal weights (α = 0.5) for VCI and TCI (blue line) against the proposed SPEI-based approach (black line) with values of α as presented in Figure 1.It may be noted that the time-series is based on values of pixels belonging to the dominant SPEI time-scale, namely 4 pixels out of 4 correspond to the 1-month SPEI scale in the case of Iberia and 8 pixels out of 16 correspond to the 9-month SPEI scale in the case of Morocco.It may be further noted that, for simplicity, we restrict to the case when VHI is computed with TCI as estimated with SMW LST.The grey vertical lines in Figure 6 identify the severe drought episodes as identified by months where the dominant time-scale of SPEI index is lower than −1.28 [66].To ensure that results are comparable, the effect of the drought episodes on the traditional and new VHI is measured based on the difference of medians normalized by the interquartile range: where p is the median of the distribution of values restricted to severe drought episodes as identified by SPEI (grey bars in Figure 6) and p , p , and p are the first quartile, the median, and third quartile of all values.This metric allows comparing the traditional VHI against the SPEI-based VHI since it assesses the capability of capturing a drought event as identified by the multi-scalar SPEI index.Results (Table 2) indicate that the SPEI-based VHI systematically presents larger absolute values of μ than the traditional VHI, indicating that the drought episodes are best represented when using the new proposed approach.To ensure that results are comparable, the effect of the drought episodes on the traditional and new VHI is measured based on the difference of medians normalized by the interquartile range: where p drought dist 50 is the median of the distribution of values restricted to severe drought episodes as identified by SPEI (grey bars in Figure 6) and p total dist 25 , p total dist 50 , and p total dist 75 are the first quartile, the median, and third quartile of all values.This metric allows comparing the traditional VHI against the SPEI-based VHI since it assesses the capability of capturing a drought event as identified by the multi-scalar SPEI index.Results (Table 2) indicate that the SPEI-based VHI systematically presents larger absolute values of µ than the traditional VHI, indicating that the drought episodes are best represented when using the new proposed approach.
Table 2. Comparison of the metric µ when using the traditional rationale to estimate VHI (µ trad ) and the SPEI-based method to estimate VHI (µ SPEI−based ) for the two regions of interest.

Discussion
Results obtained show that it is possible to estimate optimal contributions of VCI and TCI to VHI (α in Equation ( 1)) by using an independent drought index as reference.For this purpose, the multi-scalar SPEI is appropriate since it is derived using precipitation (that reflects on the moisture condition) and evapotranspiration (that reflects on the temperature condition).Based on an approach developed in a previous study [13], we show that there are indeed regions more prone to experience degraded vegetation health due to a larger contribution of either VCI or TCI.In the case of the Euro-Mediterranean window, results from the already mentioned study [13] pointed out that the contribution of moisture conditions tend to be larger over drier regions where water is a limiting factor, whereas the opposite occurs over humid regions where temperature conditions dominate, and solar energy is the limiting factor to vegetation growth.It is worth emphasizing that in Ref. [13], the approach relies on PU LST, i.e., a merged LST derived from HIRS satellite data and NCEP reanalyzed LST to fill unobserved pixels [31].However, in a climate service related to drought monitoring, other remote sensing LST products are required, as the former is not expected to be regularly updated (not at least in an operational mode).Certainly, the PU LST CDR has the advantage, when compared to other infra-red LST products (such as the two CM-SAF LST products) to be available under all-weather conditions, however the use of a unitary emissivity may also introduce a bias in the analysis.
Here, the analysis is extended by using MFG/MSG clear-sky retrievals of LST (i.e., the PMW and SMW datasets) over an area encompassing the Meteosat disk (Figure 1) and comparing the results with those obtained with PU LST.The use of these datasets to estimate the contributions of moisture (VCI) and temperature (TCI) conditions to the overall vegetation health may be an added benefit, since its practical applicability as a drought monitoring satellite product may be achieved by calibrating the contributions with the CM-SAF LST CDRs and then disseminating the product in near-real time with the currently operational MSG and future MTG.However, it is important to stress that the PU LST and CM-SAF LSTs have two fundamental differences: (1) the surface emissivity and (2) the availability under clear and all-weather conditions, respectively.The effect of emissivity may be attenuated by the fact that LSTs are used in an index such as TCI, which is a measure of temperature stress where the LST is positioned within the maximum range for a given pixel and month.Nevertheless, both effects may have an impact on the results.
The three obtained sets of parameter α (one with each LST CDR) present very similar spatial patterns (Figure 1, top panel).However, some differences are apparent like the ones between α estimated with PU LST and the two CM-SAF LSTs over the region of Algeria/Tunisia, which may be linked with the humid climate conditions found in northern Algeria [67] that may influence the retrieval of CM-SAF LSTs; or the extended patch of very low values of α over Kenya/Tanzania when using SMW LST.On the other hand, the similarity of patterns between the three datasets is worth pointing out.For instance, the east/west differences over South Africa are very well identified in all three datasets.An analysis of obtained patterns of α in terms of the different biomes (Figure 5, top panels) reveals that there are large differences between the distributions of α among the biomes.South Africa is characterized by grasslands over the western part and shrublands and forests over the eastern part.Previous studies [68,69] showed a strong relationship between NDVI and drought conditions over grasslands located in the USA and China, while Ref. [70] shows higher correlations between NDVI and SPEI over grassland regions in south Africa.Another work [71] correlates NDVI and SPEI over semiarid regions with a similar methodology used here and also shows very high correlations between the two variables over South African grasslands.
South African regions are also those where optimal correlations between VHI and SPEI present the largest improvement when compared with baseline values estimated using the traditional approach with α = 0.5 for all pixels (Figure 3).
The coherence between SPEI time-scales for correlations estimated with the three datasets is also worth noticing (Figure 1, bottom).The maps show that vegetation health responds predominantly to short drought time-scales of 1-and 3-months, which is expectable according to Ref. [70] and following the rationale that vegetation typically has a shorter memory to temperature variations.However, to the best of our knowledge, there are no works focusing on the relationship between SPEI and temperature alone.Future work on this specific topic may include a study that compares temperature (e.g., TCI) and SPEI in order to understand the drought time-scales associated to temperature conditions of vegetation.Another topic that may be tackled is the sensitivity of the contribution of α to the drought time-scale by estimating α with one time-scale for all pixels (e.g., 1-or 3-month time-scale).Although such choice will simplify the analysis, it may also introduce biases on the results, in particular since there is still no knowledge on the typical time-scales associated with temperature condition.
When comparing the three datasets, differences of about 0.2 (RMSD) were obtained between α as estimated with PU LST and CM-SAF LSTs.These differences include some pixels below (above) 0.5 when using PU LST and above (below) 0.5 when using CM-SAF LSTs.These cases represent about 19% of the total pixels.However, if cases where α jumps from 0.425 to 0.575 are dismissed, then only about 15% of the values of α estimated with PU LST and CM-SAF LSTs present relevant differences in terms of the dominant contribution.Furthermore, when comparing α estimated with the different LST CDRs there are no associated systematic errors (bias = 0).
A comparative study of severe drought episodes was finally made between the optimal VHI with SPEI-based α (VHI SPEI ) as derived from one of the LST datasets (CM-SAF SMW) and the baseline VHI estimated with α = 0.5 (VHI K ).Results obtained with the metric µ, which allows a comparison between the traditional VHI and the SPEI-based VHI by assessing the capability of capturing a drought event as identified by SPEI, point out that the SPEI-based VHI systematically presents larger absolute values and therefore it better represents the drought episodes.

Conclusions
This paper describes an approach leading to an optimal VHI for drought monitoring by identifying the contributions of VCI and TCI that maximize the correlation between VHI and a reference drought index (SPEI).The VHI and associated weight α are estimated for two clear-sky MFG/MSG LST climatological datasets (CM-SAF LSTs) and for one all-weather LST dataset based on merged HIRS observations and NCEP reanalysis (PU LST).The methodology also involves identifying the SPEI time-scales that maximize the correlation with VHI.It is shown that despite the use of datasets with different observational sources (the case of PU LST versus CM-SAF LST CDRs) or different algorithms (the case of the two CM-SAF LSTs), SPEI scales present a coherent spatial distribution over most of the pixels.Furthermore, VHI values and the respective contributions of VCI and TCI (quantified by the weight α) reveal a fair overall consistency over the studied region among datasets.As expected, the higher agreement was found between the two datasets that only differed in the LST algorithm (i.e., CM-SAF LST CDRs).Although VHI depends on both LST and NDVI, this work focuses mostly on LST since, as already pointed out, the variability of LST in terms of observation sources (e.g., HIRS versus Meteosat) and retrieval methodologies is wider than for NDVI.It is also worth recognizing that applications of NDVI in climate studies and vegetation variability assessment are currently more consolidated than applications of LST.Nevertheless, results presented here pave the grounds to the use of an optimal VHI to monitor vegetation stress and drought events from satellite data.
Indeed, results presented in this work may be viewed as a promising first step to the development of a future MSG/MTG SPEI-based VHI satellite drought monitoring based on a near-real time operational product, presenting the added benefit of taking into consideration the response in terms of vegetation health to moisture and temperature conditions over the area encompassed by the Meteosat viewing disk.Because of the expected global warming-driven increase in frequency and severity of drought episodes over regions like the Euro-and African-Mediterranean and Sahel, a Meteosat-based drought monitoring product may be of extreme importance to mitigate future socio-economic losses related to drought episodes.
Further work is needed with the aim of increasing the accuracy of the quantification of VCI and TCI share to the overall vegetation health.It would be also interesting to relate the contribution of

Figure 1 .Figure 1 .
Figure 1.Parameter α (top) estimated by iterative maximization of Vegetation Health Index (VHI) and Standardized Precipitation-Evapotranspiration Index (SPEI), and respective SPEI time-scale (bottom) for which VHI and SPEI have the largest correlation, using Thermal Condition Index (TCI) Figure 1.Parameter α (top) estimated by iterative maximization of Vegetation Health Index (VHI) and Standardized Precipitation-Evapotranspiration Index (SPEI), and respective SPEI time-scale (bottom) for which VHI and SPEI have the largest correlation, using Thermal Condition Index (TCI) computed with Princeton University Land Surface Temperature (PU LST) (left), Satellite Application Facility on Climate Monitoring (CM-SAF) Physical Mono-Window (PMW) LST (middle), and CM-SAF Statistical Mono-Window (SMW) LST (right).

Figure 2 .
Figure 2. Scatter plots (left) and histograms of differences (right) resulting from the comparison of parameter α estimated with TCI derived with CM-SAF PMW LST versus PU LST (top), CM-SAF SMW LST versus PU LST (middle), and the two CM-SAF LSTs (bottom).Bias and root mean squared deviations (RMSD) are shown on the histograms.

Figure 2 .
Figure 2. Scatter plots (left) and histograms of differences (right) resulting from the comparison of parameter α estimated with TCI derived with CM-SAF PMW LST versus PU LST (top), CM-SAF SMW LST versus PU LST (middle), and the two CM-SAF LSTs (bottom).Bias and root mean squared deviations (RMSD) are shown on the histograms.

Figure 3 .
Figure 3. Maps (top) and histograms (bottom) showing the increase in VHI-SPEI correlations, when α is changed from 0.5 to the optimum values estimated with the iterative procedure: VHI estimated with TCI derived from PU LST (left), CM-SAF PMW LST (middle), and CM-SAF SMW LST (right).

Figure 3 .
Figure 3. Maps (top) and histograms (bottom) showing the increase in VHI-SPEI correlations, when α is changed from 0.5 to the optimum values estimated with the iterative procedure: VHI estimated with TCI derived from PU LST (left), CM-SAF PMW LST (middle), and CM-SAF SMW LST (right).

Figure 4 .
Figure 4. Advanced Along-Track Scanning Radiometer (AATSR) LST biome classification version 2 (ALB-2) derived from the GlobCover classification for pixels with significant correlation between VHI and SPEI.

Figure 4 .
Figure 4. Advanced Along-Track Scanning Radiometer (AATSR) LST biome classification version 2 (ALB-2) derived from the GlobCover classification for pixels with significant correlation between VHI and SPEI.

Figure 5 .
Figure 5. Distribution of α (top panel), Δρ (middle panel), and ΔVHI (bottom panel) as estimated with PU LST (black), CM-SAF PMW LST (blue), and CM-SAF SMW LST (red) for the different biomes (identified by the labels described in Figure 4).Box and whiskers represent percentiles 10, 25, 50, 75, and 90.For each biome, the numbers indicate the respective number of pixels of the sample.

Figure 5 .
Figure 5. Distribution of α (top panel), ∆ρ (middle panel), and ∆VHI (bottom panel) as estimated with PU LST (black), CM-SAF PMW LST (blue), and CM-SAF SMW LST (red) for the different biomes (identified by the labels described in Figure 4).Box and whiskers represent percentiles 10, 25, 50, 75, and 90.For each biome, the numbers indicate the respective number of pixels of the sample.

Figure 6 .
Figure 6.Comparison between time series of VHI estimated with: the traditional method (α = 0.5) previously introduced by Kogan (blue line); and the SPEI-based method (black line) for the Iberian (top panel) and Moroccan (bottom panel) regions.Vertical grey lines in the time-series represent severe drought events as estimated with the dominant time-scale of SPEI in each region: 1-month in Iberia and 9-month in Morocco.

Figure 6 .
Figure 6.Comparison between time series of VHI estimated with: the traditional method (α = 0.5) previously introduced by Kogan (blue line); and the SPEI-based method (black line) for the Iberian (top panel) and Moroccan (bottom panel) regions.Vertical grey lines in the time-series represent severe drought events as estimated with the dominant time-scale of SPEI in each region: 1-month in Iberia and 9-month in Morocco.

Figure 7
Figure7shows the spatial distribution of VHI computed using the SPEI-based approach (VHI SPEI ) with α as in Figure1(right panels) when deriving TCI with the SMW LST.The respective spatial distribution of baseline values estimated using the approach originally proposed by Kogan (VHI K ) are also shown (left panels).The spatial patterns of the three cases clearly reflect that the SPEI-based VHI amplifies the drought episode in the two considered regions.

Figure 7
Figure7shows the spatial distribution of VHI computed using the SPEI-based approach (VHI ) with α as in Figure1(right panels) when deriving TCI with the SMW LST.The respective spatial distribution of baseline values estimated using the approach originally proposed by Kogan (VHI ) are also shown (left panels).The spatial patterns of the three cases clearly reflect that the SPEIbased VHI amplifies the drought episode in the two considered regions.

Table 2 .
Comparison of the metric μ when using the traditional rationale to estimate VHI (μ ) and the SPEI-based method to estimate VHI (μ ) for the two regions of interest.