Using NDII patterns to constrain semi-distributed rainfall-runoff models in tropical nested catchments

. A parsimonious semi-distributed rainfall-runoff model has been developed for flow prediction. In distribution, attention is paid to both timing of runoff and heterogeneity of moisture storage capacities within sub-catchments. This model 10 is based on the lumped FLEXL model structure, which has proven its value in a wide range of catchments. To test the value of distribution, the gauged Upper Ping catchment in Thailand has been divided into 32 sub-catchments, which can be grouped into 5 gauged sub-catchments where internal performance is evaluated. To test the effect of timing, firstly excess rainfall was calculated for each sub-catchment, using the model structure of FLEXL. The excess rainfall was then routed to its outlet using the lag time from storm to peak flow ( TlagF ) and the lag time of recharge from the root zone to the 15 groundwater ( TlagS ), as a function of catchment size. Subsequently, the Muskingum equation was used to route sub-catchment runoff to the downstream sub-catchment, with the delay time parameter of the Muskingum equation being a function of channel length. Other model parameters of this semi-distributed FLEX-SD model were kept the same as in the calibrated FLEXL model of the entire Upper Ping basin, controlled by station P.1 located at the centre of Chiang Mai Province. The outcome of FLEX-SD was compared to: 1) observations at the internal stations; 2) the calibrated FLEXL 20 model; and 3) the semi-distributed URBS model - another established semi-distributed rainfall-runoff NDII annual average NDII values to construct 2 alternative models: FLEX-SD-NDII Max-Min and FLEX-SD-NDII Avg ,

be very good, even better than the correlation with the NDII, because NDII does not provide good estimates during wet periods. The SWI, which is partly model-based, was not used for calibration, but appeared to be an appropriate index for validation.

Introduction
Runoff is one of the most important components of the hydrological cycle and can be monitored by the installation of a 5 gauging station. Unfortunately, there are only a limited number of high quality gauging stations available due to topographic, financial and human resources limitations. A wide variety of rainfall-runoff models have been developed in gauged and ungauged catchments in different parts of the globe. Most rainfall-runoff models are categorised as lumped models, which can provide runoff estimates only at the site of calibration. These models include FLEXL, FLEX-Topo (Euser et al., 2015;Gao et al., 2014), NAM (Bao et al., 2011;Tingsanchali and Gautam, 2000;Vaitiekuniene, 2005;Yew Gan et al., 1997), SCS 10 (Hawkins, 1990;Lewis et al., 2000;Mishra et al., 2005;Suresh Babu and Mishra, 2011;Yahya et al., 2010), and many others.
To alleviate the limitation of lumped-rainfall-runoff models, URBS was developed as a semi-distributed nonlinear rainfall runoff routing model, which can account for the spatial and temporal variation in rainfall by separating a catchment into a 15 series of sub-catchments (Mapiam and Sriwongsitanon, 2009). Therefore, URBS claims to provide runoff estimates not only at a gauging station but also at any required upstream location (Carroll, 2004;Malone, 1999). URBS has been applied successfully for real time flood forecasting in a range of catchments from small to very large basins in Australia and in many countries worldwide (Malone, 2006;Malone et al., 2003;Mapiam and Sriwongsitanon, 2009;Mapiam et al., 2014;Rodriguez et al., 2005;Sriwongsitanon, 2010). However, this model only addresses the distribution of travel times and does 20 not address the effect of distributed storage capacities that affect the partitioning of moisture and hence the water balance. Sriwongsitanon et al. (2016) proposed to use the NDII as a proxy for root zone soil moisture and showed its effectiveness in 8 sub-catchments of the Upper Ping river basins in Thailand. This is in agreement with the study carried out by Castelli et al. comparison with NDII (correlation (r) ranging from 0.951 and 0.713 with an average of 0.883). This is with the exception of some basins such as the Amazon, Murray-Darling and Mississippi (r ranging from 0.552 and 0.677) which have high percentage of forest areas, trees intercepting deep groundwater (e.g. Eucalyptus) or plenty of precipitation with low moisture stress where NDII may not correctly reflect RZWS dynamics.

5
Regarding the use of the Soil Water Index (SWI), which is partly model-based, as a proxy for root zone soil moisture, Paulik et al. (2014) found reasonable correlations between in-situ soil moisture data from 664 stations -available through the international Soil Moisture Network (ISMN) -and the SWI produced from ASCAT SSM estimates. The average of Pearson correlation coefficients was shown to be 0.54, with 64.4% of all time series greater than 0.5. SWI may be used as another index for the soil moisture state of a basin with or without moisture stress. 10 Among the wide range of existing lumped-rainfall-runoff models, FLEXL has proven to be an adequate model for runoff estimation in a wide range of catchments Fenicia et al., 2008;Gao et al., 2014;Kavetski and Fenicia, 2011;Tekleab et al., 2015). This model was further developed by Gharari et al. (2011) and Gao et al. (2016) to account for the spatial variability of landscape characteristics (FLEX-TOPO), useful for prediction in ungauged basins (Savenije, 2010). 15 Moreover, Sriwongsitanon et al. (2016) demonstrated that catchment-scale soil moisture content in the rootzone of vegetation computed from FLEXL is correlated with the remotely sensed Normalized Difference Infrared Index (NDII), as a proxy for the equivalent water thickness (EWT) in the root zone, especially during periods of moisture stress.
This study aims to utilize the fundamental model structure of FLEXL, include distributed time lags and channel routing as 20 used in URBS, and include distributed root zone soil moisture capacity per sub-catchment so as to create a new parsimonious semi-distributed FLEX model for flood and flow monitoring within the (ungauged) sub-catchments of the gauged Upper Ping River Basin. Distribution of time lags is expected to improve hydrograph shape, particularly the timing and shape of the peaks, which would improve best-fit parameters, but it does not affect the partitioning of the hydrological fluxes or the water balance. Since the root zone storage is the main control on flux partitioning, the distribution of the root zone moisture storage 25 capacity would potentially have a larger impact on model performance. Therefore, the spatial variation of the NDII, as an indicator of root zone moisture stress, has been used to distribute moisture storage capacities among sub-catchments, while the model-based SWI, as an estimator of moisture storage, was used for validation.
The main steps undertaken in the following sections are the following: 30 1. To introduce the effect of runoff timing in a catchment with multiple sub-catchments, the travel times to the outfall of each individual sub-catchment are computed on the basis of topographical indicators and the routing of the discharge from the sub-catchment outfall to stations further downstream are computed using the Muskingum method. These time lags are then applied both in the FLEX-SD model system and in the well-established URBS model, for the purpose of comparison. These two semi-distributed models only account for timing, but not for the distribution of the moisture storage capacity, a crucial parameter in runoff generation.
2. Subsequently, the effect of distribution of the root zone moisture storage is studied in the FLEX-SD model, making use of the spatial distribution pattern of the maximum and minimum range of NDII values, and the annual average NDII values to construct 2 alternative models: FLEX-SD-NDIIMax-Min and FLEX-SD-NDIIAvg, respectively. 5 3. Finally, as a validation of the model and to check if the models are capable of representing the internal moisture states, the simulated root zone moisture storage is compared to the independent data set of the Soil Wetness Index (SWI). The basin is dominated by well-forested, steep mountains in a generally north-south alignment (Sriwongsitanon and Taesombat, 2011). The areal average annual rainfall and runoff of the basin from 2001-2016 are 1,224 mm/yr and 235 mm/yr, respectively. The land use for the UPRB in 2013 can be classified into 6 main classes comprising forest, irrigated 15 agriculture, rainfed agriculture, bare land, water body, and others, which cover approximately 77.40%, 3.11%, 12.54%, 1.99%, 1.23%, and 3.73% of the catchment area, respectively (Land Development Department, LDD). The landform of the UPRB varies from an undulating to a rolling terrain with steep hills at elevations of 1500-2000 m, and valleys of 330-500m (Mapiam and Sriwongsitanon, 2009;Sriwongsitanon, 2010). Chiang Dao district, north of Chiang Mai is the origin of the Ping River, which flows downstream to the south to become the inflow of the Bhumibol Dama large dam with an active 20 storage capacity of about 9.7 billion m 3 (Sriwongsitanon, 2010). The climate of the basin is dominated by tropical monsoons.
The southwest monsoon causes a rainy season between May and October and the northeast monsoon brings dry weather and low temperatures between November and April. Only 6,142 km 2 of the total area controlled by the runoff station P.1 (situated at the centre of Chiang Mai) is selected for this study (Fig. 1). The catchment area of the station P.1 is divided into 32 sub-catchments ( Fig. 1) where the semi-distributed rainfall-runoff models are tested. 25

Rainfall data
Daily rainfall data from 48 non-automatic rain-gauge stations located within the UPRB and its surroundings from [2001][2002][2003][2004][2005][2006][2007][2008][2009][2010][2011][2012][2013][2014][2015][2016] were used in this study. These data are owned and operated by the Thai Meteorological Department and the Royal Irrigation Department. These data have been validated for their accuracy on monthly basis using double mass curve and some inaccurate data were removed from the time series before spatially averaging using an inverse distance square (IDS) to be 30 applied as the forcing data of URBS, FLEXL, and FLEX-SD. Mean areal rainfall depth for each of 32 sub-catchments varies https://doi.org /10.5194/hess-2021-598 Preprint. Discussion started: 20 December 2021 c Author(s) 2021. CC BY 4.0 License. between 1,100 (S17) and 1,402 (S11) mm/yr as shown in Figure 1 (b) while the average rainfall depth of P.1 is approximately 1,224 mm/yr.

Runoff data
The Royal Irrigation Department (RID) operates 7 daily runoff stations in the study area between 2001 and 2016 as shown in Fig. 1. Catchment P.56A was rejected from the study because it is located upstream of Mae Ngat reservoir. Outflow data 5 from the reservoir were used as input data in model calibration. Runoff data at the remaining 6 stations were used for the study since they are not affected by large reservoirs. The data have been checked for their accuracy by comparing them with average rainfall data covering their catchment areas at the same periods. Table 1 presents the catchment characteristics and hydrological data for these 6 gauging stations in the UPRB. In this study, the catchments of these 6 stations were divided into 32 sub-catchments (see Fig. 1) with areas ranging from 57 to 230 km 2 . High variation of catchment size is due to the 10 proximity between the locations of these runoff stations and the outlets of the tributaries. Runoff data have been checked for their accuracy by comparing the annual runoff coefficient between all stations. The comparison revealed that the runoff coefficients at P.20 in 2006 and 2011 are overestimated, while the runoff coefficient at P.21 in 2004 is underestimated and in 2007 and 2009 are overestimated due to incorrect rating curves (see Fig. 2). These inaccurate data would affect the results of model calibration. 15

NDII Data
The Normalized Difference Infrared Index (NDII) is a ratio of the near-infrared (NIR) and shortwave infrared (SWIR) bands, centred at 859 and 1,640 nm, respectively, as shown in Eq. (1). In this study, the NDII was calculated using 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. Atmospheric correction has been carried out to improve the accuracy and can be downloaded 20 from ftp://e4ftl01.cr. usgs.gov/MOLT (Vermote et al., 2011 (1)

SWI Data 25
The near real-time Soil Water Index (SWI) is derived from the reprocessed Surface Soil Moisture (SSM) data derived from the ASCAT sensor (Brocca et al., 2011;Paulik et al., 2014), which is a C-Band Scatterometer measuring at a frequency of 5.255 GHz in VV-polarisation (Paulik et al., 2014). The product makes use of a two-layer water balance model to describe the time series relationship between surface and profile soil moisture. This dataset of moisture conditions is available on a daily basis for eight characteristic time windows 1, 5, 10,15,20,40,60  FLEXL is a lumped hydrological model comprising five reservoirs: a snow reservoir (Sw), an interception reservoir (Si), an unsaturated soil reservoir (Su), a fast-response reservoir (Sf), and a slow-response reservoir (Ss) (Gao et al., 2014). Excess rainfall from a snow reservoir, an interception reservoir, and an unsaturated soil reservoir is divided and routed into a fastresponse reservoir and a slow-response reservoir using two lag functions. It includes the lag time from storm to peak flow (TlagF) and the lag time of recharge from the root zone to the groundwater (TlagS). Each reservoir has process equations 10 that connect the fluxes entering or leaving the storage compartment to the storage in the reservoirs (so-called constitutive functions) . The water balance equations and constitutive equations for each conceptual reservoir are summarised in Fig. 3 and Table 2. The total number of model parameters is 11. Forcing data include daily average rainfall and potential evaporation derived by the Penman-Monteith equation.

Snow reservoir 15
The snow routine, not very relevant in Thailand, can play an important role in areas with snow. When there is snow cover and the temperature (Ti) is above Tt, the effective precipitation is equal to the sum of rainfall (Pi) and snowmelt (Mi). The snowmelt (Mi) is calculated by the melted water per day per degree Celsius above Tt (FDD) (Eq. 2). The snow reservoir uses the water balance equation, Eq. (3), where Swi (mm) is the storage of the snow reservoir.

Interception reservoir 20
Interception is more important in summer and autumn. The interception evaporation Eii was calculated by potential evaporation (Epi) and the storage in the interception reservoir (Sii), with a daily maximum storage capacity (Imax) (Eqs. 4, 5). The interception reservoir uses the water balance equation, Eq. (6), presented in Table 2.

Root zone reservoir
The root zone routine, which is the core of the hydrological models, determines the amount of runoff generation. In this 25 study, we applied the widely used beta function of the Xinanjiang model (Ren-Jun, 1992) to compute the runoff coefficient for each time step as a function of the relative soil moisture. In Eq. (7), Cri indicates the runoff coefficient, Sui is the storage in the root zone reservoir, Sumax is the maximum moisture holding capacity in the root zone and β is the parameter https://doi.org/10.5194/hess-2021-598 Preprint. Discussion started: 20 December 2021 c Author(s) 2021. CC BY 4.0 License.
describing the spatial process heterogeneity of the runoff threshold in the catchment. In Eq. (8), Pei indicates the effective rainfall and snowmelt into the root zone routine; Rui represents the generated flow during rainfall events. In Eq. (9) Sui, Sumax and potential evaporation (Epi) were used to determine actual evaporation from the root zone Eai; Ce indicates the fraction of Sumax above which the actual evaporation is equal to potential evaporation, here set to 0.5 as previously suggested by Savenije (1997) otherwise Eai is constrained by the water available in Sui. The unsaturated soil reservoir uses 5 the water balance equation, Eq. (10), presented in Table 2.

Fast response reservoir
In Eq. (11), Rfi indicates the flow into the fast-response routine; D is a splitter to separate recharge from preferential flow.
Equations (12) and (13) were used to describe the lag time between storm and peak flow. Rft-i+1 is the generated fast runoff in the unsaturated zone at time t -i + 1, TlagF is a parameter which represents the time lag between storm and fast runoff 10 generation, clagF(i) is the weight of the flow in i -1 days before and Rfli is the discharge into the fast-response reservoir after convolution.
A linear-response reservoir, representing a linear relationship between storage and release, was applied to conceptualize the discharge from the surface runoff reservoir, fast response reservoirs and slow-response reservoirs. In Eq. (14), Qffi is the surface runoff, with timescale Kff, active when the storage of the fast-response reservoir exceeds the threshold Sfmax. In Eq. 15 (15), Qfi represents the fast runoff; Sfi represents the storage state of the fast response reservoirs; Kf is the timescales of the fast runoff. The fast response reservoir uses the water balance equation, Eq. (16), presented in Table 2.

Slow response reservoir
In Eq. (17), Rsi indicates the recharge of the groundwater reservoir. Equations (18) and (19) were used to describe the lag time of recharge from the root zone to the groundwater. Rst-i+1 is the generated slow runoff in the groundwater zone at time t 20 -i + 1, TlagS is a parameter which represents the lag time of recharge from the root zone to the groundwater, clagS(i) is the weight of the flow in i -1 days before and Rsli is the discharge into the slow-response reservoir after convolution. In Eq. (20) , Qsi represents the slow runoff; Ssi represents the storage state of the groundwater reservoir; Ks is the timescales of the slow runoff. The slow response reservoir uses the water balance equation, Eq. (21), presented in Table 2.

URBS model 25
URBS was developed by Queensland Department of Natural Resources and Mines in 1990 based on the structures of RORB (Laurenson and Mein, 1990) and WBNM (Boyd et al., 1987). URBS is a semi-distributed rainfall-runoff model that can provide runoff estimates not only at the calibrated station but also at the outlet of every sub-catchment at any required location upstream. The calibrated catchment area needs to be divided into sub-catchments to obtain different areal rainfall and different catchment and channel travelling time.  Table 3 presents 5 main processes used in URBS comprising the calculation of the initial loss, proportional loss, excess rainfall, catchment routing and channel routing. Excess rainfall is calculated separately between pervious and impervious areas. For the pervious area, URBS assumes that there is the maximum initial loss rate (ILmax) to be reached before any rainfall becoming the effective rainfall ( ). The initial loss (ILi) can be recovered when the rainfall rate (Ri) is less than the recovering loss rate (rlr) per time interval ( t) (see Eq. (22)). 5 Excess rainfall for each time step is calculated using Eq. (23) by weighting the excess rainfall between pervious and impervious area using a ratio of the cumulative infiltration (Fi) and the maximum infiltration capacity (Fmax). The recovering rate is included by simply reducing the amount infiltrated after every time step using the reduction coefficient (k t) as shown in Eq. (24), and the pervious excess rainfall ( ) is calculated using the Eq. (25), where pr is the proportional runoff 10 coefficient. The remaining water (1-pr) will infiltrate to the root zone storage (dFi) (see Eq. (26)). Excess rainfall is then routed to the centroid of any sub-catchment using a nonlinear reservoir relationship (Si = K ). The parameter m is the catchment non-linearity and K is the catchment travel time, which can be calculated for different sub-catchment using the multiplication between the catchment lag time coefficient (β) and square root of each sub-catchment area (A) (see Eq. (27)).

15
Thereafter, the outflow at the centroid of each sub-catchment is routed along a reach downstream of each sub-catchment using the Muskingum equation ( ℎ = Kch(XIi+(1-X)Qi)). The parameter X is the Muskingum coefficient and Kch is the channel travel time, which can be calculated for different sub-catchment using the multiplication between the channel lag coefficient (α) and the reach length (L) between the closest location in the channel to the centroid and the outlet of each subcatchment (see Eq. (28)). 20

Development of the semi-distributed FLEX model
The first step in distribution is to account for the timing of floods and the rooting of flood waves as a function of topographical factors. The resulting semi-distributed FLEX-SD model therefore is expected to better represent the shape of hydrographs, although it would not affect the partitioning of fluxes or the water balance. The root zone storage capacity is a 25 strong control on partitioning, affecting both runoff generation and evaporation. Therefore, distribution of this parameter would potentially affect overall model performance more strongly than merely the timing of the peaks. Therefore, in a second step, the NDII, as a proxy for moisture storage, is used to assess the distribution of moisture storage among subcatchments.

Accounting for distributed timing and channel-routing
FLEX-SD is set-up by applying lumped models for each sub-catchment, adding up to a semi-distributed model for a downstream calibration site. Therefore, the catchment area of any gauging station needs to be divided into sub-catchments.
Runoff estimates at each sub-catchment can be simulated using the structure of the original FLEXL by calculating different excess rainfall for each sub-catchment. The excess rainfall of each sub-catchment is routed to its outlet using the lag time 5 from rainfall to surface runoff (TlagF) and the lag time of recharge from the root zone to the groundwater (TlagS). In this study, TlagF and TlagS are calculated in hours instead of days to increase model performance. The lag time is distributed among sub-catchments using the following equations.
where, Tlag is a lag time parameter for the entire catchment of a calibrated gauging station. The lag time of each subcatchment (Tlagsub) is scaled by the square root of each sub-catchment area divided by the overall catchment area (A).
Runoff estimates from an upstream sub-catchment is later routed from its outlet to the outlet of a downstream sub-catchment using the Muskingum method (Eq. (31)) before adding to the runoff estimates of the downstream sub-catchment. 15 where, α and X are the delay time parameter and the channel routing parameter for the entire catchment, respectively. The delay time parameter of each sub-catchment (Ksub) can be calculated by the multiplication between α and the main channel length of each sub-catchment as shown in Equation (32). 20

Accounting for distributed root zone storage at sub-catchment scale using the maximum and minimum values of NDII (FLEX-SD-NDIIMaxMin model)
The Normalized Difference Infrared Index (NDII) was used to estimate root zone storage capacity for each sub-catchment.
The NDII values, which are available at 8 day intervals, were found to correlate well with the 8-day average root zone moisture content (Su) simulated by FLEXL during the dry period in eight sub-catchments in the UPRB (Sriwongsitanon et 25 al., 2016). The relation between NDII and Su can be described by an exponential function of the type: ae b(NDII) +c, with c close to zero. The maximum value that Su can achieve is Sumax, the storage capacity of the root zone. The hypothesis is that the ecosystem creates sufficient storage to overcome a critical period of drought (Gao et al., 2014;Savenije and Hrachowitz, 2017). Every year has a maximum range of storage variation. If a sufficiently long NDII record is available, then the maximum of the annual ranges of the NDII should provide an estimate of the root zone storage capacity Sumax. By 30 calibrating the hydrological FLEX model to discharge observations at the gauging stations, for each gauged catchment a https://doi.org/10.5194/hess-2021-598 Preprint. Discussion started: 20 December 2021 c Author(s) 2021. CC BY 4.0 License.
Sumax value can be calibrated. This is a representative Sumax value for a particular gauging station, consisting of n subareas, indicated by Sumaxn.
By using the NDII as proxy for root zone storage, we have developed the following equation for the proxy root zone storage capacity Sumax'i for a sub-area within a river basin consisting of 31 sub-catchments: 10 Where Sumax'i is a scaled proxy for the root zone storage capacity of each sub-catchment, and b is the remaining calibration parameter, because the constant c and the factor a of the exponential function drop out.  (35) Where Sumaxn is the calibrated value of the root zone storage capacity of the gauged catchment, and Sumax'n is the area weighted proxy for the root zone storage capacity. 20

25
Instead of applying the maximum and minimum of the annual ranges of the NDII to distribute root zone storage at a subcatchment scale, we tested the annual average NDII value of each sub-area to calculate Sumax'i as presented in the following equation.
Where NDIIi represents the annual average NDII value of each sub-catchment, while ( × → ) max and ( × → ) min indicate the maximum and minimum values of exponential function produced by the annual average NDII value within 32 sub-catchments. The parameters b and R can be determined by model calibration. The parameter R is suggested to vary between 0.2 and 0.8 to force a scaled factor Sumax'i to be more than 0 and less than 1. The average NDII value is supposed 5 to reflect the maximum moisture storage capacity as well, since a high maximum value also leads to a higher average, but is much easier to calculate. However, this method requires the introduction of the additional calibration parameter R.
In addition, all semi-distributed models were calibrated and validated at these 5 stations for a fair comparison with the results of the locally calibrated FLEXL model at each internal station (presented in Annex A). The model parameters of the 15 calibrated models were determined using the MOSCEM-UA (Multi-Objective Shuffled Complex Evolution Metropolis-University of Arizona) algorithm (Vrugt et al., 2003) by finding the Pareto-optimal solutions defined by three objective functions of the Kling-Gupta Efficiencies for high flows, low flows, and the flow duration (KGEE, KGEL and KGEE), respectively. KGEE is analysed using the following equations, where ̅ is the average observed discharge, ̅ is the average simulated discharge, is the standard deviation of observed discharge, is the standard deviation of simulated discharge, 20 and r is the linear correlation between observations and simulations. KGEL can be calculated using the logarithm of flows to emphasize low flows. The Nash-Sutcliffe Efficiency (NSE) is an independent statistical indicator, which is not utilised in the objective function but merely used to summarised model performance. The model calculates at daily time steps, but this is disaggregated to hourly to take into account the time lags. The output is again aggregated to daily time steps.  Table 4 and Table 5, respectively. Model parameters of these models are presented in Table A1. Figures 4, 5 and 6 present the output of the five models for all stations compared to observations, as accumulated flows, hydrographs on logarithmic scale, and duration curves, respectively. For comparison, the results of the calibrated and validated models at each of the 6 individual stations are presented in Figures A1, A2 and A3. Of course, these calibrated models close the water balance better, but this may be due to over-fitting. Table 4 shows that FLEXL, FLEX-SD, FLEX-SD-NDIIMaxMin and FLEX-10 SD-NDIIAvg calibrated at P.1 produce similar overall accuracy with an average NSE of 0.73, 0.73, 0.72 and 0.75, respectively during the calibration period, while URBS obtained a lower NSE of 0.68. However, FLEXL acquired higher KGE values compared to other models. while the other models were used in predictive mode. The fact that in the validation mode all semi-distributed models obtain more accurate results than the lumped and calibrated FLEXL model indicates a higher predictive capacity of the semidistributed models.
Figure 4 clearly shows that the distributed models are not capable of closing the water balance in four stations except at P.1 20 and P.67 located in the main Ping. While FLEXL, which is calibrated at each individual station can mimic the pattern, this may be due to over-fitting. In P.75, the models over-estimate the observed flow. This is due to flow regulation and water withdrawals in the managed parts of the sub-catchments (there is the large Mae Ngad dam upstream of P.75). The duration curves in Figure 6 confirm this and also show that the observed lowest flow is below the modelled flow, almost throughout. This is likely due to water abstractions for urban and agricultural water supply. In contrast, the models underestimate the 25 flows in P.21 and P.20, which are intensively used catchments with rating problems. On the other hand, the flow is overestimated at P.4A, which drains a mountainous catchment with evergreen forest. The lumped models are apparently not yet capable to distinguish well between these different landscapes. A landscape-based model as suggested by Gharari et al. (2011) and Savenije (2010) could be the next step for improvement. In general, the FLEX-SD-NDII models provides lower Sumax estimates than the other models, constraining evaporation in the dry season (which provides more realistic recessions in Figure 5), but compensates for this reduction by a smaller  value, so as to limit excessive flood generation. Since these parameters jointly control Eq. (7), they can compensate for each 15 other, leading to equifinality. If one of the parameters is constrained by additional information, as is the case here using the NDII, then this is no longer possible. The performance with respect to best fit parameters may reduce in the process, but the model has gained realism and hence predictive power.
We see that the FLEXL-SD-NDII models show the highest realism (illustrated clearly in Figure 5 and Figure A2 in the 20 appendix A) but not a very good performance in the sub-catchment P.20, although still better than the other SD models. P.20 remains a difficult sub-catchment to predict due to its flow regulation and water consumption. Also we see that adding constraints to model calibration does not always improve best-fit performance, as compared to free calibration, but that realism can be improved. To further test the realism of the models, in the following section the outputs of the models are compared to observations of NDII and the global scale SWI dataset for verification. 25 Sriwongsitanon et al. (2016) suggested that NDII can be used as a proxy for soil moisture storage in hydrology. Therefore, the 8 day average NDII values were compared to the 8 day average root zone moisture storage (Sui) as calculated by FLEXL, FLEX-SD, FLEX-SD-NDIIMax-Min and FLEX-SD-NDIIAvg. Table 6 shows the coefficients of the exponential relationships and the coefficients of determination (R 2 ) together with the NSE for the wet season, and the dry season for all six stations. The same procedure was also carried out for all 31 sub-catchments and the results shown in Table 7 Table A2. The Table confirms that FLEX-SD-NDIIAvg performs slightly better than FLEX-SD-NDIIMaxMin, but this is not surprising as it has one more calibration parameter (R), which provides an additional degree of freedom.

The relationship between the average root zone soil moisture storage (Sui) and the average NDII and SWI
It is to be noted that some sub-catchments show much lower R 2 and NSE values compared to the rest, which may be the 10 result of land use/land cover. The evergreen forest probably experiences less moisture stress compared to other land use/land cover, in which situation the NDII does not relate as well to root zone soil moisture. Therefore, Figure A5 Figure A6 presents the corresponding scatter plots for six stations and it clearly shows that the correlation is 20 much better in the dry season than in the wet season. This is not surprising, as it was argued by Sriwongsitanon et al. (2016) that the relation between NDII and root zone soil moisture can only be observed by this remote sensing product when the vegetation is experiencing moisture stress. Hence correlations between root zone soil moisture and NDII are poor during the wet season. Because the FLEX-SD-NDII was constrained by the spatial variability of NDII ranges, the good correlation between Su and NDII during the dry season may not be surprising. Therefore, an additional test was done, testing the 25 modelled Su values at daily time step with the daily SWI values, for all models.  Detailed information for 31 sub-catchments is displayed in Table 9 and Table A3. The correlation does not significantly deviate among different models for all seasons. We also show the time series plots of the average NDII (Scaling), average SWI (Scaling) and the average root zone moisture storage (Su) calculated by all models for six runoff stations in the wet and the dry seasons separately in Figures A8 and A9, respectively. One should realise, however, that the SWI is partly model 5 based and that this may affect the good correspondence during the wet season. It can therefore be concluded that the NDII is a suitable parameter to constrain hydrological models during moisture recession, but that it works less well under wet conditions. The SWI, being partly model-based, is less attractive as a model constraint, but does not suffer from a drawback during wet conditions, and hence serves well as an assessment criterion, particularly during wet conditions. As a result, the NDII appears to be useful to constrain hydrological models during dry conditions and both SWI and NDII appear to be 10 useful to test model performance and to assess moisture states of river basins.

Conclusion
Most lumped rainfall-runoff models are controlled by a gauging station at the outfall on which it is calibrated. Runoff estimation at any location upstream requires indirect approaches such as model parameter transfer from gauged stations to ungauged locations, or applying relationships between model parameters and catchment characteristics to the ungauged 15 locations. By using any of these approaches, uncertainty in runoff estimation for ungauged catchments is unavoidable. A semi-distributed hydrological model could offer a better alternative. Besides taking into account lag times and flood routing (as in FLEX-SD), it has been shown that it is required to account for the spatial variation of the moisture holding capacity of the root zone. Therefore, the model was constrained by using NDII patterns as a proxy for the spatial variation of root zone moisture leading to distributed Sumax-values among sub-catchments. We concluded that the maximum of a series of annual 20 ranges (NDIIMaxMin) and annual average (NDIIAvg) of NDII values offers an effective proxy for estimating the appropriate Sumax values in the different sub-catchments. It was shown that the two FLEX-SD-NDII models significantly improved the relationship between NDII and the modelled root zone moisture storage (Sui) of all 31 sub-catchments. Moreover, the time series of the SWI correlated very well with the modelled root zone moisture storage (Sui) of all sub-basins controlled by runoff stations. 25 The model parameters provided by the semi-distributed FLEX models are more realistic compared to the original FLEXL since they are distributed according to catchment characteristics comprising catchment area, reach length, and remote sensing indices (NDII and SWI). A next step in the analysis is to account for diversity in landscape composition and related model structures among sub-catchments , which would allow for a distinction between the main rainfall-30 runoff mechanisms belonging to different landscape types. This study confirms the result of the earlier study by Sriwongsitanon el al. (2016) who concluded that NDII can be used as a proxy for catchment-scale root zone moisture deficit https://doi.org/10.5194/hess-2021-598 Preprint. Discussion started: 20 December 2021 c Author(s) 2021. CC BY 4.0 License. when plants are exposed to water stress. However, during the wet season when soil moisture is replenished as a result of rainfall, NDII values are no longer well correlated with soil moisture. However, thepartly model-based -SWI proved to be a reliable index to estimate soil moisture both under water stressed and wet conditions.    Snow