Comparing the Normalized Difference Infrared Index (NDII) with root zone storage in a lumped conceptual model

With remote sensing we can readily observe the Earth’s surface, but direct observation of the sub-surface remains a challenge. In hydrology, but also in related disciplines such as agricultural and atmospheric sciences, knowledge of the dynamics of soil moisture in the root zone of vegetation is essential, as this part of the vadose zone is the core component controlling the partitioning of water into evaporative fluxes, drainage, recharge, and runoff. In this paper, we compared the catchment-scale soil moisture content in the root zone of vegetation, computed by a lumped conceptual model, with the remotely sensed Normalized Difference Infrared Index (NDII) in the Upper Ping River basin (UPRB) in northern Thailand. The NDII is widely used to monitor the equivalent water thickness (EWT) of leaves and canopy. Satellite data from the Moderate Resolution Imaging Spectroradiometer (MODIS) were used to determine the NDII over an 8-day period, covering the study area from 2001 to 2013. The results show that NDII values decrease sharply at the end of the wet season in October and reach lowest values near the end of the dry season in March. The values then increase abruptly after rains have started, but vary in an insignificant manner from the middle to the late rainy season. This paper investigates if the NDII can be used as a proxy for moisture deficit and hence for the amount of moisture stored in the root zone of vegetation, which is a crucial component of hydrological models. During periods of moisture stress, the 8-day average NDII values were found to correlate well with the 8-day average soil moisture content (Su) simulated by the lumped conceptual hydrological rainfall–runoff model FLEX for eight sub-catchments in the Upper Ping basin. Even the deseasonalized Su and NDII (after subtracting the dominant seasonal signal) showed good correlation during periods of moisture stress. The results illustrate the potential of the NDII as a proxy for catchmentscale root zone moisture deficit and as a potentially valuable constraint for the internal dynamics of hydrological models. In dry periods, when plants are exposed to water stress, the EWT (reflecting leaf water deficit) decreases steadily, as moisture stress in the leaves is connected to moisture deficits in the root zone. Subsequently, when the soil moisture is replenished as a result of rainfall, the EWT increases without delay. Once leaf water is close to saturation – mostly during the heart of the wet season – leaf characteristics and NDII values are not well correlated. However, for both hydrological modelling and water management, the stress periods are most important, which is why this product has the potential of becoming a highly efficient model constraint, particularly in ungauged basins.

Abstract. With remote sensing we can readily observe the Earth's surface, but direct observation of the sub-surface remains a challenge. In hydrology, but also in related disciplines such as agricultural and atmospheric sciences, knowledge of the dynamics of soil moisture in the root zone of vegetation is essential, as this part of the vadose zone is the core component controlling the partitioning of water into evaporative fluxes, drainage, recharge, and runoff. In this paper, we compared the catchment-scale soil moisture content in the root zone of vegetation, computed by a lumped conceptual model, with the remotely sensed Normalized Difference Infrared Index (NDII) in the Upper Ping River basin (UPRB) in northern Thailand. The NDII is widely used to monitor the equivalent water thickness (EWT) of leaves and canopy. Satellite data from the Moderate Resolution Imaging Spectroradiometer (MODIS) were used to determine the NDII over an 8-day period, covering the study area from 2001 to 2013. The results show that NDII values decrease sharply at the end of the wet season in October and reach lowest values near the end of the dry season in March. The values then increase abruptly after rains have started, but vary in an insignificant manner from the middle to the late rainy season. This paper investigates if the NDII can be used as a proxy for moisture deficit and hence for the amount of moisture stored in the root zone of vegetation, which is a crucial component of hydrological models. During periods of moisture stress, the 8-day average NDII values were found to correlate well with the 8-day average soil moisture content (S u ) simulated by the lumped conceptual hydrological rainfall-runoff model FLEX for eight sub-catchments in the Upper Ping basin. Even the deseasonalized S u and NDII (af-ter subtracting the dominant seasonal signal) showed good correlation during periods of moisture stress. The results illustrate the potential of the NDII as a proxy for catchmentscale root zone moisture deficit and as a potentially valuable constraint for the internal dynamics of hydrological models. In dry periods, when plants are exposed to water stress, the EWT (reflecting leaf water deficit) decreases steadily, as moisture stress in the leaves is connected to moisture deficits in the root zone. Subsequently, when the soil moisture is replenished as a result of rainfall, the EWT increases without delay. Once leaf water is close to saturation -mostly during the heart of the wet season -leaf characteristics and NDII values are not well correlated. However, for both hydrological modelling and water management, the stress periods are most important, which is why this product has the potential of becoming a highly efficient model constraint, particularly in ungauged basins.

Introduction
Estimating the moisture content of the soil from remote sensing is one of the major challenges in the field of hydrology (e.g. De Jeu et al., 2008;Entekhabi et al., 2010). Soil moisture is generally seen as the key hydrological state variable determining the partitioning of fluxes (into direct runoff, recharge, and evaporation) (Liang et. al., 1994), the interaction with the atmosphere (Legates et. al., 2011), and the carbon cycle (Porporato et al., 2004). The root zone of ecosystems, being the dynamic part of the unsaturated zone, is the key part of the soil related to numerous sub-surface processes (Shukla and Mintz, 1982). Several remote sensing products have been developed especially for monitoring soil moisture (e.g. SMOS, ERS, and AMSR-E) but until now correlations between remote sensing products and observed soil moisture at different depths have been modest at best (Parajka et al., 2006;Ford et al., 2014). There are a few possible explanations. One is that it is not (yet) possible to look into the soil deep enough to observe soil moisture in the root zone of vegetation (Shi et al., 1997;Entekhabi et al., 2010); the second is that soil moisture observations at certain depths are maybe not the right indicators for the amount of moisture stored in the root zone (Mahmood and Hubbard, 2007), which is rather determined by the vegetation-dependent, spatially variable, three-dimensional distribution and density of roots.
These mainstream methods to derive soil moisture from remote sensing have concentrated on direct observation of soil moisture below the surface. The vegetation, through the vegetation water content (VWC), perturbs this picture. As a result, previous studies have tried to determine the VWC from a linear relationship with the equivalent water thickness (EWT) that is measured by the Normalized Difference Infrared Index (NDII) (e.g. Yilmaz et al., 2008). The NDII was developed by Hardisky et al. (1983) using ratios of different values of near infrared reflectance (NIR) and short wave infrared reflectance (SWIR), defined by (ρ NIR − ρ SWIR ) / (ρ NIR + ρ SWIR ), similar to the NDVI, which is defined by discrete red and near infrared. In addition to determining the water content of vegetation, the NDII can be effectively used to detect plant water stress according to the property of shortwave infrared reflectance, which is negatively related to leaf water content due to the large absorption by the leaf (e.g. Steele-Dunne et al., 2012;Friesen et al., 2012;Van Emmerik et al., 2015). Many studies have found relationships between the EWT and reflectance at the NIR and SWIR portion of the spectrum used for deriving NDII (Hardisky et al., 1983;Hunt and Rock, 1989;Gao, 1996;Ceccato et al., 2002;Fensholt and Sandholt, 2003). Yilmaz et al. (2008) found a significant linear relationship (R 2 = 0.85) between EWT and NDII. Subsequently, they tried to determine a relationship between EWT and VWC in order to be able to correct direct moisture observations from space. However, these relationships appeared to be vegetation and crop-type dependent.
Water is one of the determinant environmental variables for vegetation growth, especially in water-limited ecosystems during dry periods. From the plant physiology point of view, water absorption from the root zone is driven by osmosis. Subsequently, water transport from the roots to the leaves is driven by water potential differences, caused by diffusion of water out of stomata, called transpiration. This physiological relationship supports the correlation between root zone soil moisture content, moisture tension in the leaves, and the water content of plants.
Hence, the root zone moisture deficit is connected to the water content of the canopy/leaves, because soil moisture suction pressure and moisture content in the leaves are directly connected (Rutter and Sands, 1958). The NDII was developed to monitor leaf water content (Hardisky et al., 1983), so one would expect a direct relation between NDII and root zone moisture deficit. The deficit again is a direct function of the amount of moisture stored in the root zone.
So, if leaf water thickness and the suction pressure in the root zone are connected, then the NDII would directly reflect the moisture content of the root zone. It would only reflect the moisture content in the influence zone of roots and not beyond that. Hence, the NDII could become a powerful indicator for monitoring root zone moisture content, providing an integrated, depth-independent estimation of how much water is accessible to roots, available for vegetation. In other words, the NDII would allow us to see vegetation as a sort of natural manometer, providing us with information on how much water is available in the sub-surface for use by vegetation. It would be an integrated indicator of soil moisture in the root zone, available directly at the scale of interest.
Thus, the hypothesis is that we can monitor the moisture content in the root zone from the observed moisture state of the vegetation by means of the NDII.
In this paper, we tested whether there exists a direct and functional relationship between a remote sensing product (the NDII) and the amount of moisture stored in the root zone, as simulated by a semi-distributed conceptual hydrological model, in which the root zone moisture content is a key state variable in the short-and long-term dynamics of the rainfall-runoff signal. Because the NDII is an indicator for water stress, the index is only expected to show a strong link with the moisture content of the root zone when there is a soil moisture deficit. Without water stress occurring within the leaves, particularly during wet periods, NDII would possibly not reflect variation in root zone soil moisture content (Korres et al., 2015).
The analysis was done using data from eight sub-basins of the Upper Ping River basin (UPRB), a tropical seasonal evergreen catchment in northern Thailand. This catchment is adequate for the purpose because it has eight well-gauged sub-basins with clearly different aridity characteristics and strong seasonality, providing a good testing ground for the comparison.
The remotely sensed NDII values have been compared to the root zone storage as modelled by a semi-distributed conceptual model (semi-distributed meaning that for each subcatchment a separate conceptual model has been used). The different sub-catchments demonstrate a variety of climatic properties that allow a more rigorous test than a fully lumped model could provide. In this way, a compromise has been found between the complexity and data requirements of a fully distributed model and the simplicity of a completely lumped model. One could argue that a fully distributed conceptual model would have been a better tool to assess the spatial and temporal pattern obtained by the NDII. This is correct, but this would have required the availability of more de- tailed spatially distributed forcing data (particularly rainfall), which were not available. Moreover, if a semi-distributed lumped model, potentially less accurate than a distributed model, provides a good correlation with NDVI, then this would be a tougher text than with a fully distributed model.  (Mapiam, et al., 2014). It has an area of approx-imately 25 370 km 2 in the provinces of Chiang Mai and Lam Phun. The basin landform ranges from an undulating to a rolling terrain with steep hills at elevations of 1500-2000 m, and valleys of 330-500 m (Mapiam and Sriwongsitanon, 2009;Sriwongsitanon, 2010). The Ping River originates in the Chiang Dao district, north of Chiang Mai, and flows downstream to the south to become the inflow for the Bhumibol Dam -a large dam with an active storage capacity of about 9.7 billion m 3 (Sriwongsitanon, 2010). The climate of the region is controlled by tropical monsoons, with distinctive dry and wet seasons and free from snow and ice. The rainy season is influenced by the southwest monsoon and brings mild to heavy rainfall between May and October. Annual average rainfall and runoff of the UPRB are approx-imately 1170 and 270 mm yr −1 , respectively. Avoiding the influence of other factors, these catchments are ideal cases to concentrate on the relationship between NDII and root zone moisture content. The land cover of the UPRB is dominated by forest (Sriwongsitanon and Taesombat, 2011).

Rainfall data
Data from 65 non-automatic rain-gauge stations covering the period from 2001 to 2013 were used. A total of 42 stations are located within the UPRB while 23 stations are situated in its surroundings. These rain gauges are owned and operated by the Thai Meteorological Department and the Royal Irrigation Department. Quality control of the rainfall data was performed by comparing them to adjacent rainfall data. For each sub-basin, daily spatially averaged rainfall, by inverse distance squared, has been used as the forcing data of the hydrological model.

Runoff data
Daily runoff data from 1995 to 2011 at eight stations located in the UPRB were adequate to be used for FLEX calibration. These eight stations are operated by the Royal Irrigation Department in Thailand. The locations of these eight stations and the associated sub-basins are shown in Fig. 1. These eight stations control the runoff of the eight sub-basins on which the eight lumped conceptual models were calibrated. Runoff data at these stations are not affected by large reservoirs and have been checked for their reliability by comparing them with rainfall data covering their catchment areas at the same periods. Catchment characteristics and available data periods for model calibration of the selected eight subbasins are summarized in Table 1.

NDII data
The satellite data used for calculating the NDII is the MODIS level 3 surface reflectance product (MOD09A1), which is available at 500 m resolution in an 8-day composite of the gridded level 2 surface reflectance products. Each product pixel contains the best possible L2G observation during an 8-day period selected on the basis of high observation coverage, low view angle, absence of clouds or cloud shadow, and aerosol loading. MOD09 (MODIS Surface Reflectance) is a seven-band product, which provides an estimate of the surface spectral reflectance for each band as it would have been measured at ground level without atmospheric scattering or absorption. This product has been corrected for the effects of atmospheric gases and aerosols (Vermote et al., 2011). The available MODIS data covering the UPRB from 2001 to 2013 were downloaded from ftp://e4ftl01.cr. usgs.gov/MOLT. The HDF-EOS conversion tool was applied to extract the desired bands (bands 2 (0.841-0.876 µm) and 6 (1.628-1.652 µm)) and re-projected into Universal Transverse Mercator (zone 47N, WGS84) from the original ISIN mapping grid.

Estimating vegetation water content using near infrared and short wave infrared
Estimates of vegetation water content (the amount of water in stems and leaves) are of interest to assess the vegetation water status in agriculture and forestry and have been used for drought assessment (Cheng et al., 2006;Gao, 1996;Gao and Goetz, 1995;Ustin et al., 2004;Peñuelas et al., 1993). Evidence from physically based radiative transfer models and laboratory studies suggests that changes in water content in plant tissues have a large effect on the leaf reflectance in several regions of the 0.7-2.5 µm spectrum (Fensholt and Sandholt, 2003). Tucker (1980) suggested that the spectral interval between 1.55 and 1.75 µm (SWIR) is the most suitable region for remotely sensed leaf water content. It is well known that these wavelengths are negatively related to leaf water content due to a large absorption by leaf water (Tucker, 1980;Ceccato et al., 2002). However, variations in leaf internal structure and leaf dry matter content also influence the SWIR reflectance. Therefore, SWIR reflectance values alone are not suitable for retrieving vegetation water content. To improve the accuracy of estimating the vegetation water content, a combination of SWIR and NIR (0.7-0.9 µm) reflectance information was utilized because NIR is only affected by leaf internal structure and leaf dry matter content but not by water content. A combination of SWIR and NIR reflectance information can remove the effect of leaf internal structure and leaf dry matter content and can improve the accuracy of retrieving the vegetation water content (Ceccato et al., 2001;Yilmaz et al., 2008;Fensholt and Sandholt, 2003). On the basis of this idea, Hardisky et al. (1983) derived the NDII: where ρ 0.85 and ρ 1.65 are the reflectances at 0.85 and 1.65 µm wavelengths, respectively. NDII is a normalized index and the values theoretically vary between −1 and 1. A low NDII value and especially below zero means that reflectance from ρ 0.85 is lower than the reflectance from ρ 1.65 , which indicates canopy water stress. The 8-day NDII values, as collected from MODIS, were averaged over each sub-basin to allow comparison to the 8day average S u (root zone storage) values extracted from the FLEX model results at each of the eight runoff stations.
We did not use field observations of soil moisture. One could argue that field observations should be used to link NDII to moisture stress. However, besides not being available, it is doubtful if point observations at fixed depth would  provide a correct measure for the moisture content in the root zone. It is more likely that vegetation distributes its roots and adjusts its root density to the specific local conditions and that the root density and distribution is not homogeneous in space and depth.
3.2 The semi-distributed FLEX model FLEX (Fig. 2) is a conceptual hydrological model with an HBV-like model structure developed in a flexible modelling framework (Fenicia et al., 2011;Gao et al., 2014a, b). The model structure comprises four conceptual reservoirs: the interception reservoir S i (mm), the root zone reservoir representing the moisture storage in the root zone S u (mm), the fast response reservoir S f (mm), and the slow response reservoir S s (mm). It also includes two lag functions representing the lag time from storm to peak flow (T lagF ) and the lag time of recharge from the root zone to the groundwater (T lagS ). Besides a water balance equation, each reservoir has process equations that connect the fluxes entering or leaving the storage compartment to the storage in the reservoirs (so-called constitutive functions). Table 2 shows the 15 equations of the FLEX model, discussed below. The 11 model parameters with their distribution values are shown in Table 3, which have to be determined by model calibration. Forcing data include the elevation-corrected daily average rainfall (Gao et al., 2014a), daily average, minimum and maximum air temperature, and potential evaporation derived by the Hargreaves equation (Hargreaves and Samani, 1985).
Splitter and Lag function 661 Table 3. Parameter ranges of the FLEX model.

Interception reservoir
The interception reservoir uses the water balance equation, Eq.
(2), presented in Table 2. The interception evaporation E i (mm d −1 ) is calculated by potential evaporation E 0 (mm d −1 ) and the storage of the interception reservoir S i (mm) (Eq. 3). There is no effective rainfall P e (mm d −1 ) as long as the S i is less than its storage capacityS i,max (mm) (Eq. 4) (de Groen and Savenije, 2006).

Root zone reservoir
The moisture content in the root zone is simulated by a reservoir (Eq. 5) that partitions effective rainfall into infiltration   (2) 0.119 (7) 0.171 (4)  and runoff R (mm d −1 ), and determines the transpiration by vegetation E t (mm d −1 ). Being the key partitioning point, the root zone storage reservoir is the core of the FLEX model. For the partitioning between infiltration and runoff, we applied the widely used beta function (Eq. 6) of the Xinanjiang model (Zhao, 1992;Liang et al., 1992), developed based on the variable contribution area theory (Hewlett and Hibbert, 1967;Beven, 1979), but which can equally reflect the spatial probability distribution of runoff thresholds. The moisture storage in the root zone reservoir is represented by S u (mm). The beta function defines the runoff percentage C r (-) for each time step as a function of the relative soil moisture content (S u / S u,max ). In Eq. (6), S u,max (mm) is the root zone storage capacity and β (-) is the shape parameter describing the spatial distribution of the root zone storage capacity over the catchment. In Eq. (7), the relative soil moisture and potential evaporation are used to determine the transpiration E t (mm d −1 ); C e (-) indicates the fraction of S u,max above which the transpiration is no longer limited by soil moisture stress (E t = E 0 − E i ).

Response routine
In Eq. (8), R f (mm d −1 ) indicates the flow into the fast response routine; D (-) is a splitter to separate recharge from preferential flow. In Eq. (9), R s (mm d −1 ) indicates the flow into the groundwater reservoir. Equations (10) and (11) are used to describe the lag time between storm and peak flow. R f (t −i +1) is the generated fast runoff from the root zone at time t −i+1; T lag is a parameter which represents the time lag between storm and fast runoff generation; c(i) is the weight of the flow in i − 1 days before; and R fl (t) is the discharge into the fast response reservoir after convolution.
The linear response reservoirs, representing linear relationships between storages and releases, are applied to conceptualize the discharge from the fast runoff reservoir, and slow response reservoir. Eq. (12) presents the water balance of the fast reservoir in which Q ff (mm d −1 ) is the direct surface runoff, with timescale K ff (d), described by Eq. (13), activated when the storage of fast response reservoir exceeds the threshold S f,max (mm), and Q f (mm d −1 ) is the fast subsurface flow, with timescale K f (d), described by Eq. (14). The slow groundwater reservoir is described by Eq. (15), which generates the slow runoff Q s (mm d −1 ) with timescale K s (d), described by Eq. (16). Q m (mm d −1 ) is the total amount of runoff simulated from the three individual components, adding up Q ff , Q f , and Q s .

Model calibration
A multi-objective calibration strategy has been adopted in this study to allow for the model to effectively reproduce different aspects of the hydrological response, i.e. high flow, low flow, and the flow duration curve. The model was therefore calibrated to three Kling-Gupta (K-G) efficiencies (Gupta et al., 2009): (1) the K-G efficiency of flows (I KGE ) measures the performance of hydrograph reproduction, especially for high flows; (2) the K-G efficiency of the logarithm of flows emphasizes low flows (I KGL ); and (3) the K-G efficiency of the flow duration curve (I KGF ) represents the flow statistics.
The MOSCEM-UA (Multi-Objective Shuffled Complex Evolution Metropolis-University of Arizona) algorithm (Vrugt et al., 2003) was used as the calibration algorithm to find the Pareto-optimal solutions defined by the mentioned three objective functions. This algorithm requires three parameters including the maximum number of iterations, the number of complexes, and the number of random samples that is used to initialize each complex. To ensure fair comparison, the parameters of MOSCEM-UA were set based on the number of model parameters. Therefore, the number of complexes is equal to the number of free parameters n; the number of random samples is equal to n · n · 10; and the number of iterations was set to 30 000. The model is a widely validated model, which is only used here to derive the magnitude of the root zone moisture storage. Therefore, validation is not considered necessary, since the model is merely meant to compare calibrated values of S u with NDII.

Deseasonalization
Seasonal signals exist both in the NDII and S u time series. This can lead to spurious correlation. Therefore, we deseasonalized both signals to eliminate this strong signal (Schaefli and Gupta, 2007) and subsequently compare the deviations from the seasonal signals of both NDII and S u . Firstly, the NDII and S u were normalized between 0 and 1. Then, seasonal patterns of NDII and S u were determined as the av-  erage seasonal signals, after which they were subtracted from the normalized data.

Spatial and seasonal variation of NDII values over the UPRB
To demonstrate the spatial and seasonal behaviour of the NDII over the UPRB, the 8-day NDII values were aggregated to monthly values for 2001-2013. Figure 3 shows examples of monthly average NDII values for the UPRB in 2004, which is the year with the lowest annual average NDII value. The figure shows that NDII values are higher during the wet season (May-October) and lower during the dry season (November-April). The lower amounts of rainfall between November and April cause a continuous reduction of NDII values. On the other hand, higher amounts of rainfall between May and October result in increasing NDII values. However, NDII values appear to vary little between July and October. The average NDII values during the wet season, the dry season, and the whole year within the 13 years are presented in Table 4. The table also shows the order of the NDII values from the highest (number 1) to the lowest (number 13). It can be seen that the annual average NDII value for the whole basin is approximately 0.  Table 5 shows the monthly averaged NDII values between 2001 and 2013 and the ranking order for each of the 14 tributaries. The results suggest that the Nam Mae Taeng, Nam Mae Rim, and Upper Mae Chaem, which have higher mean annual NDII values, have a higher moisture content than other tributaries, while Nam Mae Haad, Nam Mae Li, and Ping River sections 2 and 3, with lower mean annual NDII values, have lower moisture content than other tributaries. Monthly average NDII values for these six tributaries are presented in Fig. 4. It can be seen that during the dry season, NDII values of the three tributaries with the lowest values are a lot lower than those of the three with the highest NDII values. However, NDII values for these two groups are not significantly different during the wet season. The figure also reveals that NDII values tend to continuously increase from relatively low values in March to higher values in June. The values slightly fluctuate during the wet season before sharply falling once again when the rainy season ends, and reach their minimum values in February.

FLEX model results
Calibration of FLEX was done on the eight sub-catchments that have runoff stations. The results are summarized in Table 6. The performance of the model was quite good, as demonstrated in Table 7. In Fig. 5, the flow duration curves of runoff stations P.20 and P.21 are presented as examples of model performance. Table 7 shows the average Kling-Gupta efficiencies values for I KGE , I KGL , and I KGF , which indicate the performance of high flows, low flows, and flow duration curve for the eight runoff stations. The results for the flow duration curve appear to be better than those of the high flows and especially the low flows. However, the overall results are acceptable and can be used for further analysis in this study.

Relation between NDII and root zone moisture storage (S u )
The 8-day NDII values were compared to the 8-day average root zone moisture storage values of the FLEX model. It appears that during moisture stress periods, the relationship can be well described by an exponential function for each of the eight sub-catchments. Table 8 presents the coefficients of the exponential relationships as well as the coefficients of determination (R 2 ) for annual, wet season, and dry season values for each sub-catchment. The coefficients are merely meant for illustration. They should not be seen as functional relationships yet. The corresponding scatter plots are shown in Table 5.

P.21
Observed runoff Simulated runoff -1 -1 -1 -1 Figure 5. Examples of flow duration curves and simulated hydrographs using FLEX at runoff stations P.20 and P.21. Table 6. FLEX parameters calibrated at eight runoff stations located in the UPRB.   Fig. 6. It can be clearly seen that the correlation is much better in the dry season than in the wet season. During the wet season, there may also be short periods of moisture stress, where the exponential pattern can be recognized, but no clear relation is found when the vegetation does not experience any moisture stress. Examples of deseasonalized and scaled time series of NDII and root zone storage (S u ) values for the subcatchments P.20 and P.21 are presented in Fig. 7. The scaled time series of the NDII and S u values were calculated by dividing their value by the differences between their maximum and minimum values: NDII/(NDII max -NDII min ) and S u / (S u,max − S u,min ), respectively, while the maximum and the minimum are the values within the overall considered time series. Figure 7 shows that the scaled NDII and S u values are highly correlated during the dry season, but less so during the wet season. These results confirm the potential of NDII to effectively reflect the vegetation water content, which, through the suction pressure exercised by the moisture deficit, relates to the moisture content in the root zone. During dry periods, or during dry spells in the rainy season, as soon as the leaves of the vegetation experience suction pressure, we see high values of the coefficient of determination.
If the soil moisture in the root zone is above a certain threshold value, then the leaves are not under stress. In the UPRB, this situation occurs typically during the middle and late rainy season. The NDII then does not vary significantly while the root zone moisture storage may still vary, albeit above the threshold where moisture stress occurs. This causes a lower correlation between NDII and root zone storage during wet periods. Interestingly, even during the wet season dry spells can occur. We can see in Fig. 6, that during such a dry spell, the NDII and S u again follow an exponential relationship.
We can see that the S u , derived merely from precipitation and energy, is strongly correlated to the vegetation water observed by NDII during condition of moisture stress, without time lag (Fig. 6, and Figs. S1, S2 in the Supplement). Introduction of a time lag resulted in reduction of the correlation coefficients (see the Supplement). This confirms the direct re-sponse of vegetation to soil moisture stress, which confirms that the NDII can be used as a proxy for root zone moisture content.
The deseasonalized results of dry periods in subcatchments P.20 and P.21 are shown in Fig. 7. We found these variations of deseasonalized NDII and S u to be similar in these two sub-catchments, with the coefficients of determination (R 2 ) as 0.32 and 0.18, respectively, in P.20 and P.21. More important than the coefficient of determination is the similarity between the deseasonalized patterns. For P.20, the year 2001 is almost identical, whereas the years 2004 and 2006 are dissimilar. In general, the patterns are well reproduced, especially if we take into account the implicit uncertainties of the lumped hydrological model, the uncertainties in the 8-day derived NDII, and the data of precipitation and potential evaporation used in the model. The results of other tributaries can be found in the supplementary materials.

Is vegetation a troublemaker or a good indicator
for the moisture content of the root zone?
In bare soil, remote sensors can only detect soil moisture within a few centimetres below the surface (∼ 5 cm) (Entekhabi et al., 2010). Unfortunately, for hydrological modelling, the moisture state of the bare surface is of only limited interest. What is of key interest for understanding the dynamics of hydrological systems is the variability of the moisture content of the root zone, in which the main dynamics take place. This variability determines the rainfall-runoff behaviour, the transpiration of vegetation, and the partitioning between different hydrological fluxes. However, observing the soil moisture content in the root zone is still a major challenge (Entekhabi et al., 2010). Normally, the moisture content of the surface layer is linked to the total amount of moisture in the root zone. Knowing the surface soil moisture, the root zone soil moisture can be estimated by an exponential decay filter (Albergel et al., 2008;Ford et al., 2014) or by models (Reichle, 2008). However, the surface soil moisture is only weakly related to root zone soil moisture (Mahmood and Hubbard, 2007); it only works if there is connectivity between the surface and deeper layers, and when a certain state of equilibrium has been reached (when the short-term dynamics after a rainfall event has levelled out). It is also observed that the presence of vegetation prevents the observation of soil moisture and further deteriorates the results (Jackson and Schmugge, 1991). Avoiding the influence of vegetation in observing soil moisture (e.g. by SMOS or SMAP) is seen as a challenge by some in the remote sensing community (Kerr et al., 2001;Entekhabi et al., 2010). Several algorithms have been proposed to filter out the vegetation impact (Jackson and Schmugge, 1991), also based on NDII (e.g. Yilmaz et al., 2008   In this study, we found that vegetation, rather than becoming a problem, could become key to sensing the storage dynamics of moisture in the root zone. The water content in the leaves is connected to the suction pressure in the root zone (Rutter and Sands, 1958). If the suction pressure is above a certain threshold, then this connection is direct and very sensitive. We found a highly significant correlation between NDII and S u , particularly during periods of moisture stress. During dry periods or dry spells in the rainy season, as soon as the leaves of the vegetation experience suction pressure, we see high values of the coefficient of determination. Observing the moisture content of vegetation provides us with direct information on the soil moisture state in the root zone. We also found that there is almost no lag time between S u and NDII. This illustrates the fast response of vegetation to soil moisture variation, which makes the NDII a sensitive and direct indicator for root zone moisture content. In fact, the canopy acts as a kind of manometer for the root zone moisture content.

The validity of the hypothesis
In natural catchments, it is not possible to prove a hypothesis by using a calibrated model. There are too many factors contributing to the uncertainty of results: the processes are too heterogeneous, the observations are not without error, the climatic drivers are too uncertain and heterogeneous, and finally, there is substantial model uncertainty, both in the semidistributed hydrological model and in the remote sensing model used to determine the 8-day NDII product. In this case, we have selected a lumped conceptual model, which is good at mimicking the main runoff processes, but which lacks the detail of distributed models. Distributed models, however, require detailed and spatially explicit information (which is missing) and are generally over-parameterized, turning them highly unreliable in data-scarce environments. On top of this, there is considerable doubt if they provide the right answers for the right reasons.
This paper is not a modelling study but a test of the hypothesis whether the observed NDII correlates with the modelled root zone storage. We have seen in Fig. 6 that the correlation is strong during periods of moisture stress, but that when the root zone is near saturation the correlation is weak. But we also saw that even in the wet season, during short dry spells, the correlation is strong. Even when the seasonality is removed, the patterns between NDII and S u in Fig. 7 are similar, although there are two dry seasons when this is less the case (in 2004 and 2006). So given the implicit uncertainty of the hydrological model, the uncertainty of the meteorological drivers, as well as the river discharges to which the models have been calibrated, and the uncertainty associated with the relationship between NDII and EWT, the good correspondence between the NDII and the root zone storage of the model during periods of moisture stress support the potential value of the NDII as a proxy for root zone storage in a conceptual model. It is in our view even likely that the differences between the signals of the NDII and the S u are rather related to model uncertainty, the uncertainty of the climatic drivers, the uncertainty in the relationship between NDII and EWT, and the problems of determining accurate NDII estimates over 8-day periods, than due to a weak correlation between the root zone storage and the NDII.

Implication in hydrological modelling
Simulation of root zone soil moisture is crucial in hydrological modelling (Houser et al., 1998;Western and Blöschl, 1999). Using estimates of soil moisture states could increase model performance and realism, but moreover, it would be powerful information to facilitate prediction in ungauged basins (Hrachowitz et al., 2013). However, until now, it has not been practical (e.g. Parajka et al., 2006;Entekhabi et al., 2010). Assimilating soil moisture in hydrological models, either from top-soil observation by remote sensing, or from the deeper soil column by models (Reichle, 2008), is still a chal-lenge. Several studies showed how difficult it is to assimilate soil moisture data to improve daily runoff simulation (Parajka et al., 2006;Matgen et al., 2012).
There are several reasons why we have not compared our results with soil moisture observations in the field. Firstly, observations of soil moisture are not widely available. Moreover, it is not straightforward to link classical soil moisture observations to the actual moisture available in the root zone. Most observations are conducted at fixed depths and at certain locations within a highly heterogeneous environment. Without knowing the details of the root distribution, both horizontally and vertically, it is hard, if not impossible, to estimate the water volume accessible to plants through their root systems. We should realize that it is difficult to observe root zone soil moisture even at a local scale. But measuring root zone soil moisture at a catchment scale is even more challenging. State-of-the-art remote sensing techniques can observe spatially distributed soil moisture, but what they can see is only the near-surface layers if not blocked by vegetation. The top layer moisture may or may not be correlated with the root zone storage, amongst others, depending on the vegetation type, but it is definitely not the same.
By observing the moisture content of the leaves, the NDII represents the soil moisture content of the entire root zone, which is precisely the information that hydrological models require as this is the component that controls the occurrence and magnitude of storage deficits and thereby the moisture dynamics of a system. This study clearly shows the temporal correlation between S u and NDII. From the relationship between NDII and S u , we can directly derive a proxy for the soil moisture state at the actual scale of interest, which can potentially be assimilated in hydrological models. Being such a key state variable, the NDII-derived S u could become a potentially powerful and otherwise unavailable constraint for the soil moisture component of hydrological models. This would mean a breakthrough in hydrological modelling as it would allow a robust parameterization of water partitioning into evaporative fluxes and drainage even in data-scarce environments. Given the implicit uncertainties in hydrological modelling, this new and readily available proxy could potentially enhance our implicitly uncertain modelling practice. More importantly, it would open completely new venues for modelling ungauged parts of the world and could become extremely useful for discharge prediction in ungauged basins (Hrachowitz et al., 2013).
We should, of course, be aware of regional limitations. The proxy only appears to work for periods of moisture stress. This study considered a tropical seasonal evergreen ecosystem, where periods of moisture stress regularly occur. In ecosystems which shed their leaves or go dormant, other conditions may apply. We need further investigations into the usefulness of this approach in catchments with different climates. In addition, the phenology of the ecosystem is of importance, which should be taken into consideration in follow-up research. Finally, a comparison with distributed or Hydrol. Earth Syst. Sci., 20, 3361-3377, 2016 www.hydrol-earth-syst-sci.net/20/3361/2016/ semi-distributed models would be a further test of the value of the NDII as proxy for the root zone moisture content.

Conclusions
The NDII was used to investigate drought for the UPRB from 2001 to 2013. Monthly average NDII values appear to be spatially distributed over the UPRB, in agreement with seasonal variability and landscape characteristics. NDII values appear to be lower during the dry season and higher during the wet season as a result of seasonal differences between precipitation and evaporation. The NDII appears to correlate well with the moisture content in the root zone, offering a potential proxy variable for calibration of hydrological models in ungauged basins.
To illustrate the importance of NDII as a proxy for root zone moisture content in hydrological models, we applied the FLEX model to assess the root zone soil moisture storage (S u ) of eight sub-catchments of the UPRB controlled by eight runoff stations. The results show that the 8-day average NDII values over the study sub-basin correlate well with the 8-day average S u for all sub-catchments during dry periods (average R 2 equals 0.87), and less so during wet spells (average R 2 equals 0.61). The NDII appears to be a promising proxy for root zone moisture content during dry spells when leaves are under moisture stress. The natural interaction between rainfall, soil moisture, and leaf water content can be visualized by the NDII, making it an important indicator both for hydrological modelling and drought assessment.
The potential of using the NDII to constrain model parameters (such as the power of the beta function β, recharge splitter D, and C e in the transpiration function) in ungauged basins is an important new venue, which could potentially facilitate the major question of prediction in ungauged basins.

Data availability
The data set can be found at: https://zenodo.org/record/ 60491.
The Supplement related to this article is available online at doi:10.5194/hess-20-3361-2016-supplement.