Effect of Using Multi-Year LULC and Monthly LAI Inputs on the Calibration of a Distributed Hydrologic Model

This study explores the use of satellite-based LULC (Land Use / Land Cover) data while simultaneously correcting potential evapotranspiration (PET) input with Leaf Area Index (LAI) to increase the performance of a physically distributed hydrologic model. The mesoscale hydrologic model (mHM) was selected for this purpose due to its unique features. Since LAI input informs the model about vegetation dynamics, we incorporated the LAI based PET correction option together with multi-year LULC data. The Globcover land cover data was selected for the single land cover cases, and hybrid of CORINE (coordination of information on the environment) and MODIS (Moderate Resolution Imaging Spectroradiometer) land cover datasets were chosen for the cases with multiple land cover datasets. These two datasets complement each other since MODIS has no separate forest class but more frequent (yearly) observations than CORINE. Calibration period spans from 1990 to 2006 and corresponding NSE (Nash-Sutcliffe Efficiency) values varies between 0.23 and 0.42, while the validation period spans from 2007 to 2010 and corresponding NSE values are between 0.13 and 0.39. The results revealed that the best performance is obtained when multiple land cover datasets are provided to the model and LAI data is used to correct PET, instead of default aspect-based PET correction in mHM. This study suggests that to minimize errors due to parameter uncertainties in physically distributed hydrologic models, adequate information can be supplied to the model with care taken to avoid over-parameterizing the model.


Introduction
Decision-making about the sustainability of water resources requires supporting tools to manage the resources effectively. Hydrological models are useful tools to capture spatiotemporal variations in hydrological fluxes in many basins. A model consists of various parameters that define the characteristics of the basin. Hydrologic models have become progressive tools for studying the effects of anthropogenic activities and environments on hydrology and ecology [1]. Good modeling should be characterized by a high degree of confidence in the simulated outputs and minimization of uncertainty in model results. These uncertainties may be due to data input integrity, measured output data, non-optimal parameters, and model bias.
The conventional knowledge is that the best simulation is associated with less uncertainty in its output. The best model is the one that gives results close to reality with the use of the least parameters and model complexity, i.e. a parsimonious model. Hydrologists are faced with improving the quality of model outputs by using various procedures to reduce these uncertainties. However, errors due to model bias are difficult to address and are usually mentioned by the modeler when reports are made. To address uncertain-ties relating to non-optimum parameter values, model calibration is employed. This procedure is targeted at obtaining optimum parameter values for the model setup in the specific basin, such that the error due to non-optimal parameters is insignificant compared with errors due to data integrity [2].
Various strategies have been adopted to calibrate hydrological models. Scientists have introduced series of mathematical functions that guide parameter optimization during the calibration of the model. Hydrologists are also engaged in exploring model setups and data inputs to obtain the best model performances. Various resolutions of Digital Elevation Models (DEM), Land Use Land Cover (LULC) source and resolutions, soil map resolution, and data length have been studied to obtain their effects on model performance. These responses vary among models due to model types and structures [3]. Physically distributed hydrological models are capable of simulating multi-layered hydrological processes. Rather than utilizing extensive hydrological and meteorological data during calibration, physically distributed hydrologic models require the evaluation of many parameters that describes the physical characteristics of the catchment [4]. These models provide output upon which decisions about effective water resources management are based especially with spatiotemporal changes in LULC. Therefore, the accuracy of these outputs is of high importance to hydrologist, hence the need for efficient calibration procedure to design simulations representing the interested basin.
Bouslihim et al. [5] studied the effect of input resolution on the Soil and Water Assessment Tool (SWAT) model's performance. Soil data from different databases with different resolutions were used to calibrate the model, showing that the data's quality and resolution affect the estimation of hydrological fluxes, especially soil water storage and water yield. Ficchì et al. [6]also studied the impact of temporal resolution of input data on model performance using different precipitation time-steps to simulate short duration hydrologic events like flash floods. Storm duration, flood hydrograph peaks, rainfall-runoff lag time and temporal precipitation variability, are the most influential variables for model performance at different time steps. The effect of calibration data and length on model performance was observed by Ilampooranan et al. [7], using crop yield data to calibrate a SWAT model rather than the usual streamflow data to identify a deficiency in model structure. The research highlighted the importance of additional data sources to improve the robustness and predictive ability of distributed models. Li et al. [8] further studied the effect of calibration length on lumped hydrologic model performance. They showed that longer calibration periods do not necessarily result in better model performance and optimum parameter values can be attained with fewer calibrations in humid catchments.
The use of a multi-objective function to calibrate the model to obtain improved model performance was studied by Sahraei et al. [9]. They found that although the performance of the models was improved when calibrated with multi-objective functions, they were unable to identify all aspects of simulation errors relating to the model uncertainty. This led to the use of segmentation-based model-wrapper to group the models together and significantly increased the performance compared to individual models. The influence of farming activities was also studied by Gao et al. [10], where the impacts of dominant crop rotation techniques on SWAT hydrologic model performance were observed. It was concluded that proper sets of crop rotations in the model could improve model performance within the rotation period. However, applying the same set of crop rotations for longer periods in the model decreases the model performance, and the crop rotation pattern often leads to model over-parameterization and equifinality.
Orth et al. [11] studied the influence of model complexities on performance. The study concluded that complexity does not necessarily cause increased model performance. Model performance varies differently with the hydrological variable under consideration. The transfer of information from gauged catchments to ungauged catchments for the prediction of runoff is common in hydrologic modeling, which justifies the need to observe the response of models to varying methods of achieving this. Yang et al. [12] observed the influence of regionalization methods on the performance of hydrological models in different climatic regions and showed that model parameters vary among regionalization methods and climatic regions with more variations associated with the change of climatic region, which affects the performance of the models.
Land cover depicts the physical state of the land surface, in terms of the spatial distribution of soil, water, vegetation and anthropogenic activities on the earth surface (Parveen et al. [13]; Prasad & Ramesh [14]; Risal et al. [15]. Leaf Area Index (LAI) is an essential parameter of plant community that describes the total area of leaves per unit ground area and is directly related to the amount of light available for interception by a plant. It defines the total one-sided area of leaf tissue per unit ground surface area [16] and influences transpiration, interception of rainfall and sunlight, gas exchange and photosynthetic processes [17], through the canopy structural characteristics [18].
The sensitivity of vegetation parameters to hydrological processes was assessed by Das et al. (2018) using Land Use Land Cover (LULC) changes in the eastern India basin between 1995 and 2005 to set up VIC model. The LAI emerged as the most sensitive parameter that influences hydrological fluxes. LULC changes through deforestation, urbanization, expansion of croplands led to reduction in the extent of canopy cover for interception and transpiration, which affects hydrologic performance. This necessitated the study of the impact of multiple LULC datasets and resolutions on hydrologic performance observed by Jin et al. [1] and suggested using a high-resolution LULC dataset for modeling when a year LULC dataset is available. Further, the use of 10-year LULC dataset to calibrate the SWAT model in the study improved model performance by 2.2% -13.9%.
The use of yearly datasets identifies the LULC changes better by smoothing deviations over the input data of the semi-distributed model. However, little information exists about the effect of multiple LULC datasets on the performance of fully distributed physically based hydrologic models like the mesoscale Hydrologic Model (mHM). Hence, this study aims at 1) exploring the performance of a distributed hydrologic model by calibrating the hydrologic model using a series of land use land cover datasets and 2) analyzing the influence of leaf area index, which is a significant characteristic of forest areas on the performance of the model. This research involves different model configuration cases with varying LULC, potential evapotranspiration (PET) setup, and comparing the performance and hydrograph of each simulated output.

Study Area
The study area is the Karasu basin located at the Euphrates River's headwaters in Turkey ( Figure 1). The basin is mostly dominated by pasture, grass and bare land. The total area is about 10,275km 2 . This area is hugely mountainous with an elevation range of 1125-3500m and mean basin slope of 20%. Karasu Basin boundaries are within the longitudes 38 o 58' E to 41 o 39' E and latitudes 39 o 23' N to 40 o 25'N. Studies showed that 60-70% of the total annual runoff volume comes during snowmelt season. The annual mean precipitation of the basin fluctuates from 400 to 450mm/yr [20]. Streamflow data for the Karasu basin are collected by General Directorate of Turkish State Hydraulic Works (DSI) at the Karasu Karagöze and Aşağı Kağdarıç gauging station [21]. The basin downstream is characterized by huge dams that make accurate runoff estimation crucial to efficient water resources allocation for flood control, hydropower generation, irrigation, and water supply [22].

Model Description
The mHM is a fully distributed hydrologic model based on numerical approximations of hydrologic processes such as canopy interception, snow accumulation and melting, soil moisture dynamics, infiltration and surface runoff, evaporation, subsurface storage and discharge generation, deep percolation and base flow, and discharge attenuation and flood routing [23]. The model discretization in this setup is based on grid cells of 0.0625 o spatial resolution, with each grid generating runoff that is routed using adaptive timestep, constant celerity routing method. The mHM utilizes a Multiscale Parameter Regionalization technique to account for sub-grid variability of catchment morphology (Demirel et al. [24], Dembélé et al. [25], Höllering et al. [26]) which enables the swapping between different scales while calculating all fluxes and routing streamflow on a preferred model scale. The model contains 62 global parameters that can be adjusted during the calibration process and pedo-transfer functions for soil parameterization [24]. Samaniego et al. [27] can be consulted for explicit description of model formulation and parameter description. The mHM was preferred as the interest model due to its three unique features 1) multiscale parameter regionalization technique, which links basin's physical characteristics to parameter value via pedo-transfer functions 2) ability to process multiple LULC data 3) use of leaf area index (LAI) data not only for calculating interception loss but also for correcting PET.

Data Preparation
The basic data for setting up the mHM can be classified into meteorological data, morphological data, land cover data and gauging station data. The essential meteorological variables are the precipitation and average air temperature at hourly to daily temporal resolution. However, based on a strategy of estimating potential evapotranspiration, additional data like net radiation, wind speed, the absolute vapor pressure of air, potential evapotranspiration, maximum and minimum air temperature. The morphological variables are the digital elevation model, soil maps with textural properties, geological maps containing specific yield, permeability and aquifer thickness. For this study, the E-OBS dataset was utilized with the precipitation, average air temperature data and potential evapotranspiration obtained from Rakovec et al. [28] with a resolution of 0.25 o and daily time steps, input directly into the model, hence no need for its estimation in the model ( Table 1). All morphological data are provided at a resolution of 0.001953 o and the meteorological data at 0.25 o . The model is run at a daily time-step with and outputs are internally resampled to a spatial resolution of 0.0625 o by the model.

Land Cover Data
The Coordination of Information on the Environment (CORINE), MODIS and GLOBCOVER land cover data were used. The Coordination of Information on the Environment (CORINE) Land Cover (CLC) is the first land cover map in Europe based on photographic interpretation of satellite images for the years 1990, 2000, 2006, 2012 and 2018, producing the CLC 1990, CLC 2000, CLC 2006, CLC 2012 and CLC 2018 products respectively. The spatial resolutions of the products are 100m and 250m, [29] with the land cover inventoried into 44 classes [30]. The main advantage of the CORINE Land Cover (CLC) inventory is its frequent update by European countries, although the level of detail of the source data is a limitation [30]. This limitation suggests that the CLC inventory is a useful resource for land cover change studies on the regional scale. The MODIS Land Cover collection is a global land cover data set designed to support research related to current, seasonal to the decadal status of global land properties, with 17 land cover classes at a 500m spatial resolution and yearly temporal coverage with a 71.6% overall accuracy [29]. The GLOBCOVER dataset is a global project initiated by European Space  [31]. The project entailed the generation of world land cover map using 300m Medium Resolution Imaging Spectrometer (MERIS) instrument aboard ENVISAT [32]. The mHM identifies only three land cover classification. This classification includes forest (coniferous and deciduous), impervious (urban and built-up areas, water bodies, and consolidated soils), and pervious (grasslands, croplands, and bare soil) [23]. Hence, the land cover data were reclassified to suit the data requirements.

Sensitivity Analysis
Hydrologic models are associated with uncertainties, which have always prompted researchers better to understand the hydrologic process and better representation in models. These models are associated with some level of complexities that vary among models to represent necessary variables correctly. Sensitivity Analysis (SA) is a method employed by modelers to deal with these complexities during calibration by analyzing the contribution of each model parameter to the uncertainty in model prediction. The basis for SA is an effect called the "Sparsity of Factors" principle. This principle, which originates from the field of Statistical Design of Experiments, states that when a process is influenced by several variables, its behavior is likely driven by a small subset of these variables [33]. These subsets of variables characterize sources of uncertainties in model prediction and upon which research towards model performance improvement should be focused [34]. SA can summarily be applied for assessment of model similarities [35], obtaining regions of sensitivity [36], factor and model reduction [37] and uncertainty apportionment [38]. In this study, we used PEST Tool [39] to assess most important parameters based on Jacobian matrix to include in the model calibration.

Model Calibration and Validation
In this study, four different cases (Table 2) based on input data were created to achieve the set objectives. These cases contain different configurations and representations of hydrologic processes. The Case 1 (C1) model was configured using one-year GLOBCOVER land cover data and no correction of PET data using LAI information in the region. The PET data was only subjected to correction with the aspect of the grid cells. The Case 2 (C2) model was also configured with the same land cover data but with LAI correction. This strategy allows the PET data to be corrected based on the LAI data supplied to the model. The Case 3 (C3) model was configured with multiple land cover data from a different database. It contained 1990 and 2000 CORINE land cover data and consecutive MODIS land cover data from 2001 to 2008 for the calibration period. The PET data, in this case, is also subjected to only the aspect correction and ignores the LAI correction. The Case 4 (C4) model contains the same land cover configuration as C3; however, the PET data were subjected to LAI correction. The model was calibrated with the discharge data from the gauging station from 1990 to 2006 while 2007 to 2010 was used for the validation.

Dynamically Dimensioned Search Algorithm (DDS)
The DDS was developed to obtain best global solutions within a model evaluation limit based on heuristic global search algorithm [40]. The DDS algorithm begins its search globally and becomes localized as the iteration numbers approaches the limit of function evaluations. This is achieved by decreasing the number of dimensions in the neighborhood dynamically and probalistically, and further distribute the reducing probability equally to each decision variable [41]. The need to modify parameters simultaneously during calibration to preserve the current gain in calibration result was the motivation behind the algorithm. The critical feature of DDS was triggered by the poor result associated with manual calibration of watersheds especially with progress in the calibration process. It becomes necessary to calibrate one or few parameters simultaneously as calibration result improves in order to prevent loss of previously gained result [40].

Model Performance
The Nash-Sutcliffe efficiency (NSE) was used to evaluate the integrity of the model simulations for each scenario. This objective function quantifies the accuracy of prediction by comparing the differences between predicted and observed discharge values. NSE is a scaled metric that compares the expected variable with the observed variable's mean as the naïve forecast. An efficiency of zero, indicates that using the observed data as forecast would give the same error as the predicted data, a negative efficiency means that the mean of the observed data is a better predictor than the model prediction [42]. The inverse (NSEinv) and logarithm (NSEln) of this metric was also adopted to capture the flow dynamics. Pushpalatha et al. [43] reported that NSEln and NSEinv are useful metrics for the evaluation of low flows, although the latter shows no sensitivity to high flows unlike the NSEln, which shows little sensitivity, and NSE, which is highly sensitivity to high flows. The Kling-Gupta Efficiency (KGE) was also adopted to evaluate the annual flows [44]. KGE consists of three components targeted at addressing the flaws observed in NSE. This evaluation criterion is based on NSE decomposition into the linear correlation between observed and simulated flows, the flow variability term and biased term [45]. The optimum value is 1 and unlike NSE, KGE does not have an inherent benchmark against which flows are compared and the value to be compared is subjective to the modeler [46]. Other evaluation criteria used are PBIAS and Mean Absolute Error (MAE), which measured the simulated values' average tendency to be larger or smaller than their observed ones and the absolute difference between observed and simulated flows respectively. The optimum value of MAE is 0 and low magnitude values indicates accurate model prediction. Positive values of PBIAS indicate some overestimation of predicted data while negative values indicates model underestimation of predicted values.. [47] suggested the use of some hydrological indices for the evaluation of model simulation. For this study, the long-term annual average discharge, flow in dry season, flow in wet season and base flow were selected as indices to evaluate the modeled cases. The wet season was taken from May to October while the rest of the months were considered as dry season. The base flow was estimated as 30% of the average daily flow and the performance was evaluated using the procedure alighted in Chirachawala et al. [48].

Land Use Land Cover Variations
LULC maps in physically distributed hydrological models are essential inputs for efficient modeling of spatial variations in hydrologic fluxes, excess rainfall and infiltration in particular. This study adopted the CORINE and MODIS LULC dataset and the changes in the land cover characteristics shown in Figure 2 using the mHM recognized land cover classes. The variation in land cover classes between the upper and lowest land cover map adopted is also shown in Figure A1. Using the conventional CORINE land cover nomenclature to characterize the land cover in the basin, the most dominant land cover obtained in the study period is the natural grasslands, which increased from 718.63km 2 in 1990 to 929.63km 2 in 2012, although a slight decline was noticed in 2000 before the increase in subsequent years. A significant increase in the area covered by water bodies was also observed. In 1990, 0.04% of the watershed, which is 1.1km 2 , is covered by water. This increased to 11.17km 2 in 2000 and varied annually until it reached 9.42km 2 in 2012. Road networks in the basin occupied the least land cover in the studied maps covering 0.013% of the total basin area in 1990 and expanded to 0.054% in 2012. Artificial surfaces like buildings, roads, industry and land cover comprising of urban fabrics increased by 7.35% between 1990 and 2012. This may be responsible for the decrease in the agricultural areas from 1194.53km 2 to 1139km 2 , indicating about 4.7% decline. In the same vein, forest areas declined by 6.3% between 1990 and 2012. Glacier and perpetual snow in the watershed were observed to cover more area in 2012. In 1990, the snow area was around 434.62km 2 and steadily increased annually in the observed maps to 473.98km 2 in 2012. However, the hydrologic model of interest was designed to recognize three land cover classifications; forest, impervious and pervious. MODIS has no separate forest class in the study domain after reclassification but more frequent (yearly) observations than CORINE, hence these two land cover products complement each other. Sensitivity analysis in this study was used to identify model parameters sensitive to discharge prediction in the basin using the Parameter Estimation and Uncertainty Analysis (PEST) tool [39]. PEST performs all calibration tasks, including sensitivity and uncertainty analysis, upon the user needs. The tool allowed the identification of the sensitive parameters to discharge prediction in the model set up and these were normalized using the maximum sensitivity value. Of the 62 global parameters in the model, 20 parameters were selected as most sensitive as shown in Table A1 and they were afterwards used to calibrate the hydrologic model.

Sensitivity Analysis
The mHM parameters are divided into groups based on the processes they influence. The soil moisture parameters were observed to be very sensitive to the model output, although this sensitivity varies among the parameters. Eleven of these parameters out of about twenty soil moisture parameters were sensitive in the setup. This actually depends on the method used for modeling soil water movement, the Feddes equation was used in this study. The sensitivity of PET parameters differs with choice of PET estimation method. For this study, PET correction using LAI and aspect maps were eventually adopted. The PET parameters related to correction using LAI parameters were found to be sensitive. The process of interflow was also found to be sensitive to streamflow predictions in the model set up. This was shown in the sensitivity analysis characterized with 3 out of the 5 interflow parameters showing huge sensitivity. Lastly, the percolation process, through the recharge coefficient parameter was found to be sensitive to the streamflow dynamics in the study area. These sensitive parameters were selected to calibrate the model to obtain optimum parameter values with 3000 iterations.

Model Calibration
All four cases were calibrated with the sensitive parameters obtained from the sensitivity analysis. These cases are used to consider the simultaneous effect of PET correction technique and multiple LULC maps. The estimated NSE value ranges between 0.226 and 0.422, with the highest in Case 4 and the lowest at Case 1 respectively. PBIAS and MAE value for Case 1 and Case 3 were the highest, indicating a better performance, with Case 2 and Case 4 the better cases. KGE values for the modeled cases during calibration ranges between 0.419 and 0.544 for the calibration, with Case 4 having the highest KGE value indicating the best performance with this metric. NSEln and NSEinv were adopted to filter out the high flow sensitivity of conventional NSE during estimation. The range of NSEln was -0.36 and 0.05 while that of NSEinv was-0.65 and -0.61. Considering the individual effect of PET correction method on the model performance, Case 1 and Case 2 were configured with the same input data but different PET correction method. Case 3 and Case 4 were also configured with the same input data but different PET correction method. In this vein, the NSE value of Case 2 revealed higher value than that of Case 1 while that of Case 4 also revealed a higher NSE value than that of Case 3. To isolate the effect of multiple land cover data on the model, Case 1 and Case 3 were configured with the same input dataset but different land cover configuration. Case 2 and Case 4 were also configured with the same input data but different land cover configuration. In this comparison, Case 3 and Case 4 proved better with higher NSE values. KGE values for these cases also follows the same pattern with the Case 4 having the best performance among the cases as shown in Table 3.

Model Validation
For the validation phase, the NSE value obtained were slightly lower than the calibration phase ranging from 0.19 and 0.39, with the Case 4 also showing the best performance. The NSEln metric revealed its highest value as Case 2 and the lowest with Case 3, recording efficiency of 0.26 and -0.07 respectively. KGE also recorded its highest value in Case 2 and lowest in Case 4. However, Case 4 performed better under MAE, with the lowest value among the observed Cases as shown in Table 3. The hydrograph of the observed and simulated average monthly discharges for all the cases is shown in Figure 3 and the cumulative distribution frequency of daily discharge data is shown in Figure A2. It should be noted that all cases underestimated the frequency of the discharge until it reached about 1000m 3 /s, when Case 2 and Case 4 overestimated the frequency, as did the other cases when the discharge accumulated to about 1400m 3 /s.

Hydrological Indices
The annual discharge, wet season flows, dry season flows and base flow of the modeled cases were used to evaluate the model's performance by comparing them with the observed data from the gauging station, as shown in Table 4. For annual discharge, wet season flows and dry season flows, Case 2 performed better than the others, although Case 4 also performs optimally. For annual discharge, the percentage difference of Case 2 is 2.89% while that of the worst performance is Case 3, with a percentage difference of 44.37%. The dry season flow performance showed a different trend, with Case 3 performing better than the other cases, with a percentage difference of 1.62% against the lowest performing case which was Case 2 with a difference of 22.7%.
To show the effect of misclassification of land cover on discharge, monthly discharge in the basin was predicted by providing the model with CORINE 2006 and MODIS 2006 sequentially. The CORINE land cover product characterized with forest, pervious and impervious land cover predicted higher discharges between May and October, while the MODIS land cover product with just pervious and impervious classes simulated higher discharges for the rest of the months.

Discussion
The need to estimate the impacts of climate and land use change on discharge regime has made hydrologic models popular worldwide, especially with their ability to predict flows at gauges and ungauged catchments. However, these predictions are subjected to uncertainty due to model bias, input data errors, and errors in model parameter values [49]. The goal of hydrologic modelers is to reduce uncertainties associated with the model outputs upon which decision about hydrologic fluxes are based. There have been numerous studies focusing on different ways to improve model performance by exploring various input data characteristics that affects the predictions in conceptual and physically distributed models [7]- [11], [50]. Models like mHM are spatially distributed since they consist of equations involving one or more space coordinates, which are used for simulation of spatial variation in hydrological variables within a catchment as well as simple outflows and bulk storage volumes. Such models make considerable demands in terms of computational time and data requirements and are costly to develop and operate [2]. The results presented in this study provide the spatial model's response to varying characteristics of the input data based on the objective under study. Recall that the scenarios labeled as Case 1 to Case 4 and their varying input characteristics. Case 1 and Case 2 were configured with single LULC data, with LAI influence on the Case 2 only. The same applies to Case 3 and Case 4 but with multiple land cover classes, which makes it possible to observe LAI's effect on the model performance. LAI is a key biophysical vegetation property for interception that describes biome-specific canopy structure and an essential variable in evaluating the interrelation of the components of the water-soil-plant-atmosphere system, which influences environmental studies [51]. The vegetation of an area due to its LAI is a significant component of land surface models that influence evapotranspiration simulation and, consequently, affects base flow, recharge, infiltration excess, saturation excess, subsurface storm flow and catchment wetness [52]. As hypothesized, the cases LAI data performed better than their counterparts during calibration and validation process. The implication is that predictions are more accurate when leaf area index is included in the model for correction on evapotranspiration. Precise evaluation of evapotranspiration is a vital aspect of water balance estimation. About 60% of terrestrial precipitation returns to the atmosphere by plant transpiration (40%), or through direct soil evaporation (20%) [53]. The correct simulations of LAI seasonal dynamics and stomatal aperture in an eco-hydrological model are prerequisite of good simulation of canopy radiation exchanges and transpiration fluxes which are important components of the accurate water balance estimation [54].
Furthermore, the impact of single and multiple land cover data was observed. The result associated better model performance with the cases with multiple land cover data. This corroborates Jin et al. [1], where the use of multiple land-use datasets improved model performance although with some complexities in simulation complexity due to greater number of LULC patches. LULC maps are essential input of physically distributed models since they spatially delineate the basin's different morphological characteristics that decide water's fate on the earth's surface. A detailed characterization of land cover is essential for urban areas, as coarse spatial description might result in biased estimates by hiding urban heterogeneous evapotranspiration [55]. Petrucci & Bonhomme [56] also found a similar result to our study where increasing geographical information clearly improves model performances with the land use classification providing the highest benefit. Finally, considering the combined effect of both morphological characteristics, the case with multiple land cover data and LAI to correct evapotranspiration performed better than the other cases and the case with single land cover data, without LAI correction performed worse. It is noteworthy that the discharge data were used in the model's calibration, which leaves the possibility of variation in model performance if different calibration method or data was used. The data-intensive nature of physically distributed models often leads to over-parameterization while effectively representing each hydrologic process that constitutes such a model. This opens up another phase of research in hydrologic modeling to solve over-parameterization and identify strategies to improve model performance without overwhelming the model with information.

Conclusions
Hydrologic modeling remains an important decision support tool for efficiently allocating water resources and solving global water security issues, especially with the climate's non-stationarity and changing land use. The model outputs are characterized with some degree of uncertainty by underestimating or overestimating the simulated values. The best models are thus characterized with reduced uncertainties, and different methods are being adopted to achieve this. This study explored the impact of LULC maps and LAI on the performance of a fully distributed hydrologic model. Various model cases were set up to evaluate the independent effects of LAI and LULC map on model performance. The result showed that using LAI for accurate estimation of evapotranspiration values adopted for the simulation of hydrologic processes increased the model performance when the model was calibrated using discharge data. In addition to this, providing mHM with multiple LULC datasets has better performance than those with single LUCL data. This study reiterated the observation of Jin et al. [1] about the impact of multiple land cover datasets on model performance. We further suggest using optimum land cover dataset when available to enable the model to capture recent and relevant LULC patches that will be utilized in the simulation of the hydrologic processes. Especially in rapidly developing landscapes, this will be necessary to capture the basin dynamics. Furthermore, in the absence of multiple or yearly LULC maps, other strategies that will reduce model uncertainty like using LAI data for accurate estimation of evapotranspiration can be adopted.  Table 1. Karasu basin model set-up data were kindly provided by Oldrich Rakovec from UFZ in Leipzig, Germany. Discharge data is publicly available at the repository of State Hydraulic Work (DSİ) https://www.dsi.gov.tr/Sayfa/Detay/744 and daily data from the gauging station was compiled and then provided by Dr Muhammet Yılmaz. Due to data license restrictions, the observed geographical, hydrological and meteorological data cannot be redistributed from the authors, but can be obtained from Rakovec et al. [28], respective organizations and/or UFZ in Leipzig, Germany upon request. The model calibration results can be retrieved from the corresponding author upon request. The source code of the mHM is publicly available at https://doi.org/10.5281/zenodo.4575390.