Spatio-temporal fusion of NDVI data for simulating soil water content in heterogeneous Mediterranean areas

ABSTRACT Recent studies have demonstrated that the soil water content (SWC) of Mediterranean ecosystems can be simulated by combining ground data and remote sensing observations of Normalized Difference Vegetation Index (NDVI). The application of this approach in heterogeneous and fragmented areas, however, requires the use of spatio-temporal fusion (STF) methods to properly account for the actual NDVI variability of the examined ecosystems. One of these methods, which was specifically developed to produce annual NDVI data series in Mediterranean regions, is currently applied to MODIS and TM/ETM+ images taken over a highly fragmented green urban area in Florence (Central Italy). The performances of this STF method, called SEVIS, are indirectly evaluated by comparing local SWC measurements to simulations driven by the original (MODIS) and fused (MODIS+TM/ETM+) NDVI datasets. The results obtained confirm the critical dependence of the applied SWC simulation strategy on the efficient accounting for the actual NDVI evolution of the observed ecosystem. In particular, the use of the fused NDVI dataset corrects almost completely for the strong SWC underestimation produced by the original MODIS images during the summer dry period, significantly improving all accuracy statistics (r 2 from 0.564 to 0.855, RMSE from 0.101 to 0.044 cm3 cm−3 and MBE from −0.046 to 0.000 cm3 cm−3).


Introduction
Mediterranean areas are characterized by the presence of prolonged summer dryness which can cause stress to both semi-natural and agricultural vegetation. The analysis of long-term temperature and precipitation data series highlights that water stress periods occur more frequently and with higher intensity with respect to the past (Giorgi & Lionello, 2008), enforcing the idea that Mediterranean areas are among the most sensitive to climate change (Scarascia-Mugnozza, Oswald, Piussi, & Radoglou, 2000). This increases the importance of predicting and monitoring water availability and fluxes through semi-natural and agricultural ecosystems in sight of promoting sustainable water uses for different needs.
Actual evapotranspiration (ET A ), defined as the sum of plant transpiration and soil evaporation, is a major term of site water budget which is highly variable in space and time depending on numerous environmental factors (Olivera-Guerra, Mattar, & Galleguillos, 2014). The estimation of daily ET A on regional to local scales is therefore a critical issue which has been addressed by the use of several techniques (Wilson, Hanson, Mulholland, Baldocchi, & Wullschleger, 2001). Remote sensing data, in particular, have been widely applied to this scope, thanks to their capability to provide precious information on important environmental features (e.g. soil and vegetation type and condition), at various spatial and temporal scales. Two approaches are generally applied for the estimation of ET A through remote sensing techniques, i.e. those based on energy balance and water balance methods (e.g. Glenn, Huete, Nagler, Hirschboeck, & Brown, 2007;Pereira, Allen, Smith, & Raes, 2015). Both approaches present pros and cons, but most operational applications have been performed through water balance methods driven by remotely sensed vegetation indices (VI), among which the normalized difference vegetation index (NDVI) is the most widespread (Glenn, Nagler, & Huete, 2010). This is the case for the method proposed by Maselli, Chiesi, Angeli, Papale, and Seufert (2014), named NDVI-Cws, which was specifically aimed at estimating ET A in water limited Mediterranean areas. The method developed is based on the combination of meteorological and satellite NDVI data and can be applied to drive a simplified model which simulates the daily variability of soil water content (SWC) (Battista et al., 2018).
The operational application of this simulation approach, however, is often hampered by the high spatial heterogeneity and fragmentation which are typical of many agricultural and semi-natural ecosystems. This situation is common in many regions, requiring NDVI datasets with improved spatial and temporal features which can properly describe the highly changeable phenological stages of local vegetation. While such datasets are currently acquired by new satellite systems (i.e. Sentinel 2), their availability for retrospective analyses depends on the use of spatio-temporal fusion (STF) methods, which can combine different image sets to yield synthetic products having the desired spatial and temporal features (Emelyanova, McVicar, van Niel, Li, & van Dijk, 2013;Gao, Masek, Schwaller, & Hall, 2006). Numerous STF methods have therefore been proposed in the last few decades, reaching various levels of efficiency and applicability (see Chen et al., 2015, for a review). Most of these methods estimate reflectance or VI changes relying on the definition of similarities between high spatial resolution (HR) pixels (e.g. Gao et al., 2006), which is suboptimal to yield long-term projections in areas where the main vegetation types show diverging phenologies.
To address this issue, Maselli, Chiesi, and Pieri (2016) proposed a new method based on the estimation of pixel similarities through the decomposition of low spatial resolution (LR) NDVI multitemporal data series. The method, refined and definitively named SEVIS (Spatial Enhancer of VI image Series) in Maselli, Chiesi, and Pieri (2018) is specifically suited to predict the diversified NDVI evolutions of ecosystems characterized by high spatio-temporal variability and complexity. SEVIS could therefore find its maximum utility to guide the estimation of ET A and the simulation of SWC in fragmented and heterogeneous Mediterranean areas. The current paper aims at ascertain this hypothesis concerning an urban park in Florence (Central Italy) having spatially fragmented and phenologically diversified vegetation cover. The performances of the modelling approach driven by different NDVI datasets are evaluated against SWC measurements taken during a study year (2009).

Study area and data
The selected study site is located within the urban Park of Cascine, which is the largest green area in the city of Florence (Central Italy) (Figure 1). The park lies along the Arno river, at an altitude of about 40 m a.s.l. The climate is sub-humid, with mean annual rainfall of about 800 mm and mean annual temperature of 15.7°C; rainfall is concentrated in autumn and spring, while summer is usually hot and dry. The park is covered by broadleaved and coniferous trees, grasses and some buildings interspersed with small agricultural areas. With the exception of few large meadows, the areas grown with grasses and annual crops (i.e. mainly vegetables) are highly fragmented and surrounded by vine-yards and olive groves. The soil is sandy-clay-loam (52% sand, 24% silt and 24% clay) and has a depth over 2 m at the study site.
A fully equipped agrometeorological station was placed in a field of about 0.15 ha placed near the centre of the park (43.7854°N, 11.2183°E) and usually grown with non-irrigated grasses (Figure 1 (b)). The station collected standard daily meteorological data (i.e. minimum and maximum air temperature, precipitation and solar radiation) from July 2008 to February 2013. The station also collected daily SWC measurements by means of a capacitive probe (Theta Probe ML2X) placed at 0.2 m depth.
LR MODIS NDVI images, which are collected by sensors on the Terra and Aqua satellites, were freely downloaded from the USGS database (http://lpdaac. (MVC) images have a 250-m pixel size and are composited over 16-day periods. Two cloud-free Landsat scenes of 200 × 200 pixels were utilized as HR data. The first scene was taken by Landsat 7 ETM+ on 10/04/2009; this scene was not significantly affected by the ETM+ scan line corrector failure due to its proximity to the central part of the frame. The second scene was taken by Landsat 5 TM on 29/07/2009 (Figure 1(b)). Both scenes had a spatial resolution of 30 m and were downloaded from the US Geological Survey website in a pre-processed (i.e. georeferenced) format.

SWC simulation strategy
The methodology applied for the simulation of SWC at field scale is schematized in Figure 2 and briefly described in the next paragraphs. The methodology uses ancillary information, daily meteorological datasets (i.e. air temperature, precipitation and solar radiation) and NDVI images which must be descriptive of the actual vegetation conditions and relevant seasonal evolutions in the study site. More details are provided in the mentioned references together with full justifications of all choices and assumptions made.
The simulation of the water budget is based on a simplified one-layer bucket model which considers only the main inputs and outputs to the soil. Through this model, the volumetric SWC (VSWC) of day i is calculated as: where P i is precipitation, O i is outflow (i.e. percolation plus runoff), and ET Ai is evapotranspiration, all referred to day i. The model requires the definition of soil field capacity and saturation, which are jointly used to bound the upper VSWC variability; the lower VSWC boundary is usually set to zero. A full description of this modelling approach with relevant advantages and limitations is given in Battista et al. (2018). The daily ET A estimates required by Equation (1) are obtained through the NDVI-Cws method (Maselli et al., 2014). This method uses NDVI values, converted into fractional vegetation cover, to account for the contribution of site transpiration and evaporation by means of the following equation: where ET 0 is potential evapotranspiration, FVC is the fractional vegetation cover, and Cws and AW are the two short-term water stress scalars, all referred to day i.; Kc Veg and Kc Soil correspond to the maximum crop coefficient for vegetation and soil, respectively. The two water stress scalars are predicted through the combination of precipitation and ET 0 , while the other terms are derived from the literature (see Maselli et al., 2014, for details).

STF of NDVI data
In retrospective cases, NDVI images having HR and temporal frequency can be obtained only by the application of STF methods such as SEVIS. As previously mentioned, this method had been specifically conceived to cope with the diverging phenologies of spatially heterogeneous Mediterranean vegetation types based on few single-date HR images. Thus, SEVIS should be particularly suitable for the estimation of the annual NDVI evolutions in small urban areas where artificial and semi-natural land cover types are mixed. SEVIS consists of the following five steps: (i) application of the Sequential Maximum Angle Convex Cone (SMACC) algorithm (Gruninger, Ratkowski, & Hoke, 2004) to the LR multitemporal dataset for identifying NDVI endmembers of vegetation types having uniform phenological evolutions and relevant abundance images; (ii) use of this information to drive a fuzzy maximum likelihood (ML) classification of the available HR imagery; (iii) spatial enhancement of the SMACC abundance images based on the HR ML membership grade images; (iv) modification of the LR multitemporal NDVI endmembers to cope with the intra-class variability of vegetation features which can be found in relatively large areas; (v) recomposition of the multitemporal NDVI endmembers and the spatially enhanced abundance images to produce synthetic HR imagery. When applied to 250-m 16-day MODIS MVC and 30-m Landsat TM/ETM+ images, SEVIS produces 23 spatially enhanced NDVI images for each year.

Data processing
The described methodology was applied to simulate the water budget of the study site during the year 2009 using the available satellite and ancillary data. First, the MODIS NDVI multitemporal dataset and the two Landsat scenes were pre-processed as described in Maselli, Papale, Puletti, Chirici, and Corona (2009). An annual NDVI profile was then directly extracted from the 250-m MODIS pixel corresponding to the study site. Next, SEVIS was applied to the MODIS multitemporal images and to the two Landsat scenes, producing 23 spatially enhanced images from which a new annual NDVI profile was derived for the same site. Both LR and HR NDVI profiles were temporally interpolated to daily values, linearly converted into FVC and used to feed Equations (1) and (2) (see Battista et al., 2018, for details).
The meteorological variables needed by these equations were derived from the records of the agrometeorological ground station; few missing data (30 days) were replaced by data interpolated at 250m spatial resolution, which were not significantly different from the measurements. ET 0 was estimated by applying the method of Jensen and Haise (1963). The crop coefficients Kc Veg and Kc Soil were set to 1.2 and 0.2, respectively (Chiesi et al., 2013).
When applying Equation (2), outflow was equalled to the water exceeding the maximum soil water holding capacity corresponding to a rooting depth of 0.6 m, a field capacity of 0.33 and a saturation of 0.44. The first value was estimated taking into account the existing ecosystem type, while the second two were obtained by applying a pedotransfer function to soil texture (Jabloun & Sahli, 2006).
The daily SWC estimates obtained using both the LR (MODIS) and HR (synthetic) NDVI profiles were finally evaluated versus the available measurements, summarizing the results by means of common accuracy statistics (i.e. the coefficient of determination, r 2 , the root mean squared error, RMSE, and the mean bias error, MBE).

Results
The 2009 meteorology of the study site is typical of a Mediterranean climate, with moderate rainfall total (760 mm) and marked summer dryness (Figure 3). ET 0 , which is maximum just after the summer solstice, is much higher than rainfall, reaching a total of about 1140 mm. This gives rise to a water stress period from July to September, which, however, is only partly reflected in the NDVI profile extracted from the MODIS pixel corresponding to the ground station ( Figure 4).
These NDVI values, in fact, are relatively high and nearly stable during the whole 2009, showing only minor variation between 0.4 and 0.6. This obviously corresponds to a nearly stable FVC and is reflected into the estimated ET A , which peaks in June-July similarly to ET 0 and has a high annual total (around 500 mm). The likely ET A overestimation is confirmed by the marked SWC underestimation which is visible in Figure 5 during the dry summer period.
Light on this issue is cast by the visual examination of the Landsat images, which indicate the mixed nature of the MODIS pixel corresponding to the ground station (Figure 1). Such observation is confirmed by the results of the SMACC algorithm, which identifies five main multitemporal NDVI endmembers roughly corresponding to mixed forests, urban and industrial lands, residential areas, irrigated orchards and gardens, and grasslands ( Figure 6). The respective LR abundance images indicate that all these cover types, with the exception of urban areas, are present in the MODIS pixel of the ground station. More particularly grasslands, which is the actual vegetation class surrounding the station, represents only 15% of this pixel, and therefore only marginally contributes to its NDVI signal.
The applied STF method modifies the LR abundance images by introducing relevant HR detail. Two examples of synthetic NDVI images produced for mid-April (spring maximum) and mid-August (summer minimum) are shown in Figure 7. The spatial heterogeneity of the park, particularly around the ground station, leads to extreme NDVI variability especially in summer.
The spatially enhanced abundance images entirely attribute the HR pixel of the station to grassland. Consequently, the NDVI evolution derived from the synthetic images is similar to that of endmember 5 and typical for this cover type (Figure 4): NDVI ranges from a spring peak  of 0.6 to a summer minimum of 0.3. The notable summer NDVI reduction implies a corresponding decrease in estimated ET A of about 100 mm, which leads to an annual total around 400 mm. The effect on SWC is again shown in Figure 5: both visual analysis and accuracy statistics indicate a notable increase in the accordance between SWC measurements and estimates (r 2 improves from 0.564 to 0.855, RMSE from 0.101 to 0.044 cm 3 cm −3 and MBE from −0.046 to 0.000 cm 3 cm −3 ). These findings are confirmed by a similar accuracy assessment carried out considering separately the four seasons. The results reported in Table 1 indicate that the improvement in SWC simulation is negligible in the first part of the year, while it is particularly evident during the summer dry  Figure 7. Synthetic NDVI images of mid-April (a) and mid-August (b) produced by the applied STF method. The white line and the red square are the same as in Figure 1; the NDVI ranges from 0 (black) to 1 (white). period; the improvement is prolonged in the last part of the year.

Discussion and conclusions
All methods currently applied have been put forward and tested in previous publications, where a full discussion of relevant performances and limitations is provided. The present investigation concerns the combination of these methods aimed at a specific purpose, i.e. the simulation of site water budget and SWC in a spatially fragmented area within the main urban park in Florence. The simulation of daily SWC requires ET A estimates descriptive of the actual vegetation status in the examined site. The currently applied NDVI-Cws ET A estimation method is a refinement of classical algorithms which, based on the Kc concept, use vegetation indices to correct ET 0 . The method is robust and easily applicable to a variety of ecosystem types, but critically depends on the capacity of the used NDVI data to describe the quantity and functions of local vegetation. This creates problems for the retrospective application of the method in spatially fragmented and heterogeneous Mediterranean areas, where the main vegetation types show multiple growing cycles which cannot be fully characterized by any single NDVI dataset (Maselli et al., 2016). In these cases, the inaccuracy which is introduced in the estimation of ET A can deteriorate the simulation of site water balance and SWC, as was actually found in the current experiment.
Such situation calls for the use of STF methods which can efficiently integrate the spatial and temporal features of different satellite datasets over mediumlong periods. Several STF methods have been proposed for the estimation of NDVI changes (see Chen et al., 2015, for a review), but they generally provide suboptimal performances for long-term projections in areas characterized by diverging phenological evolutions . The currently applied STF method, SEVIS, specifically addresses this issue, being capable of producing annual synthetic NDVI image series with improved spatial detail based on a singledate HR image. An exhaustive inter-comparison of this method with other, more classical algorithms was performed in the same investigation . The results of this test indicated that SEVIS outperforms the other algorithms in the presence of fragmented vegetation types showing diversified growing cycles (i.e. winter and summer crops, evergreen and deciduous forests, etc.). This is particularly the case for long-term projections, i.e. when predicting the NDVI values of seasons temporally distant from the acquisition date of the HR image.
The complex nature of the currently investigated urban area is confirmed by the asynchronous seasonal evolutions of the vegetation classes identified by decomposing the LR multitemporal NDVI dataset. Contemporaneously, the visual analysis of the used HR imagery indicates the extremely heterogeneous spatial patterns of these classes, which compose a complex mosaic of artificial and semi-natural land surfaces. In these cases, SEVIS generally finds its maximum utility, characterizing the seasonal evolution of fragmented Mediterranean vegetation types more efficiently than the original LR imagery. This expectation is actually supported by the grassland NDVI values currently obtained from the MODIS and synthetic images; while the former are strongly influenced by the presence of other vegetation types (annual and perennial crops, forest, etc.), the latter are clearly representative of nonirrigated Mediterranean herbaceous vegetation. In particular, the synthetic NDVI values better reproduce the expected annual development of the study grassland in response to the main driving meteorological factors, i.e. the spring peak due to optimal growing conditions and the summer decline following water shortage.
Such differences obviously produce different predictions of grassland ET A during the study year. The NDVI-Cws method driven by the synthetic images simulates a higher spring peak and much lower summer values, which can be considered as a more realistic consequence of the mentioned climatic factors. This is finally reflected into an enhanced SWC simulation capacity, which is quantitatively demonstrated versus relevant daily measurements. In particular, the improvement is mostly due to a better SWC simulation during the dry summer period, when the difference between the MODIS and synthetic NDVI evolutions is particularly evident.
Similar results are expectable in all other situations characterized by spatially fragmented vegetation cover, particularly in Mediterranean areas having asynchronous seasonal NDVI evolutions. Further investigations should be conducted to confirm this expectation in other regions having different land cover and vegetation patterns.