Using Remotely Sensed Information to Improve Vegetation Parameterization in a Semi-Distributed Hydrological Model (SMART) for Upland Catchments in Australia

Appropriate representation of the vegetation dynamics is crucial in hydrological modelling. To improve an existing limited vegetation parameterization in a semi-distributed hydrologic model, called the Soil Moisture and Runoff simulation Toolkit (SMART), this study proposed a simple method to incorporate daily leaf area index (LAI) dynamics into the model using mean monthly LAI climatology and mean rainfall. The LAI-rainfall sensitivity is governed by a parameter that is optimized by maximizing the Pearson correlation coefficient (R) between the estimated and satellite-derived LAI time series. As a result, the LAI-rainfall sensitivity is smallest for forest, shrub, and woodland regions across Australia, and increases for grasslands and croplands. The impact of the proposed method on catchment-scale simulations of soil moisture (SM), evapotranspiration (ET) and discharge (Q) in SMART was examined across six eco-hydrologically contrasted upland catchments in Australia. Results showed that the proposed method produces almost identical results compared to simulations by the satellite-derived LAI time series. In addition, the simulation results were considerably improved in nutrient/light limited catchments compared to the cases with the default vegetation parameterization. The results showed promise, with possibilities of extension to other hydrologic models that need similar specifications for inbuilt vegetation dynamics.


Introduction
Terrestrial vegetation plays an important role in the water and energy cycles. Changes in the vegetation cover affect land surface properties, carbon generation and consumption, and significantly impact regional and global climate systems dynamics [1][2][3]. For the hydrological cycle, vegetation density (e.g., leaf area index (LAI)) and physiological properties (e.g., stomatal conductance) control partitioning of rainfall into runoff and evapotranspiration [4]. Specifically, stomatal conductance controls leaf-scale transpiration rates, and LAI controls transpiration rates and canopy interception losses. As a result, changes in LAI not only impact evapotranspiration rates, but also can affect subsequent hydrologic processes including soil moisture, baseflow and runoff [5,6]. Such long term changes in vegetation cover in response to environmental factors can also manifest themselves through nonstationary hydrologic responses where long-term trends in annual runoff to precipitation ratios in 20 anthropogenically unaffected catchments in Australia are observed [7].
Given the importance of vegetation in the hydrologic cycle, appropriate representation of the vegetation processes (i.e., photosynthesis, respiration, carbon allocation, and phenology) in hydrologic models is desirable. Physically based ecohydrologic models such as the Regional Hydro-Ecological Simulation System (RHESSys) [8] and Tethys-Chloris [9] have been developed to simulate ecohydrological processes at catchment scales. These ecohydrologic models integrate hydrological and ecological processes of a catchment where vegetation dynamics are internally simulated as a function of water and energy availability, and vegetation is used as a state variable for solving water/energy balances [10,11]. As the application of ecohydrologic models is limited by a large number of parameters and high computational demand, it is still a challenge to effectively consider dynamic vegetation effects on the water and energy balances in hydrological models [12]. Some physically-based hydrologic models assume vegetation as a static component of the hydrologic cycle or prescribe the same monthly LAI climatology (climatologically representative value of LAI by taking a long-term average of monthly LAI values) year after year to incorporate vegetation processes [13]. For example, in a distributed hydrology-vegetation model, Wigmosta et al. [14] defined maximum canopy interception storage capacities and resistance based on seasonal LAI values. In a semi-distributed variable infiltration capacity (VIC) model, Liang et al. [15] parameterized land cover types based on monthly LAI climatology values, canopy resistance and relative fractions of roots in the soil layers. In the HYDRUS software package that simulates water, heat, and solute transport in porous media [16], potential evapotranspiration is partitioned to transpiration and soil evaporation according to the Beer's law that attenuates light through the canopy based on light extinction coefficient and LAI [17].
These limited parameterization approaches using constant or inter-annually invariant monthly LAI values have been shown to impact accuracy of simulated hydrologic fluxes. Ford and Quiring [18] showed that incorporating inter-annual LAI dynamics from the moderate resolution imaging spectroradiometer (MODIS) in the VIC model improved soil moisture simulations particularly during dry periods. However, the improvement has not been consistent. Donohue et al. [12] emphasized explicit representation of vegetation dynamics in the Budyko model particularly when it is applied at small timescales and/or small catchments as vegetation dynamics affect annual and seasonal vegetation water use. For this reason, many efforts have been made to include vegetation component in the Budyko model framework [19][20][21][22]. Wegehenkel [23] coupled a hydrologic model with a static and a dynamic vegetation module implemented at a catchment scale, and showed that simulated evapotranspiration and groundwater recharge are significantly different from those derived from temporally invariant vegetation conditions. Tang et al. [24] also showed that incorporating observed MODIS LAI compared to LAI climatology has larger impacts on evapotranspiration compared to soil moisture particularly in regions with larger LAI inter-annual variability. Tesemma et al. [4] also illustrated that VIC simulated runoff from the Goulburn-Broken catchment in Australia is improved when observed monthly LAI data from the Global Land Surface Satellite (GLASS) are used instead of long-term mean monthly LAI climatology.
As incorporating LAI dynamics at large scale depends on the quality of remotely sensed observations and data availability, empirical approaches have been developed to include LAI dynamics in hydrologic models while keeping model structure unchanged. To incorporate changes in LAI dynamics for climate change impact assessment, Tesemma et al. [25] examined predictive capability of a set of three-parameter empirical functions that relate annual or monthly LAI to a range of climate variables (precipitation, potential evapotranspiration, soil moisture and temperature) for three different land covers types in the Goulburn-Broken catchment in Australia. The best non-linear empirical LAI model calibrated against remotely sensed LAI data were first applied to simulate changes in LAI under climate change scenarios. Results of future runoff simulations using VIC model illustrated significant differences in the mean annual runoff when estimated LAI from the empirical function were utilized instead of LAI climatology [26]. Mohaideen and Varija [27] also used the VIC model for simulating various hydrologic fluxes for a catchment that has experienced significant land cover changes over a 10-year period. Satellite-derived LAI and land cover data were used for grid-wise vegetation parameterization, and results showed significant increases/decreases in hydrologic fluxes for different land cover types.
To incorporate vegetation dynamics in hydrologic models without explicitly simulating carbon cycle, this study aims to improve the LAI parameterization method in a semi-distributed hydrologic model, called the Soil Moisture and Runoff simulation Toolkit (SMART) [28]. SMART is a semi-distributed hydrologic modeling toolbox for computationally efficient hydrologic simulations. The computational efficiency is achieved by formulating two-dimensional representative hillslopes or equivalent cross sections (ECSs) in every sub-basin where each ECS consists of topographically defined hydrologic response units (HRUs). These delineations will result in a highly efficient semi-distributed hydrological modeling framework compared to distributed hydrologic models that perform model simulations at very grid cell [29,30]. Despite SMART's computational efficiency, its vegetation parameterization needs to be improved. The current version of SMART relies on a simple method proposed by Vaze et al. [31], where monthly LAI climatology is proportionally scaled by the ratio of monthly rainfall to monthly rainfall climatology for three land cover classes (pasture, crop, and tree). However, the accuracy of this method on simulating hydrologic fluxes and states has not been fully assessed compared to prescribing inter-annual variability of LAI from remotely sensed data. Therefore, we aim to develop a viable alternative to the existing LAI parameterization in SMART using remotely sensed LAI data to improve catchment-scale simulations of soil moisture (SM), evapotranspiration (ET) and streamflow (Q). We also expect that the proposed framework can be conveniently incorporated in any hydrologic model where simple parameterization of vegetation dynamics is desired. The proposed LAI scaling method was evaluated by simulating soil moisture, evapotranspiration, and streamflow over six ecohydrologically contrasted catchments located in high elevation mountain ranges in Australia. For this research objective, this paper is organized as follows. In Section 2, we describe the study area and data used in this study together with a summary of the exiting LAI model in SMART, and the proposed LAI parameterization method. In Section 3, we present and discuss the evaluation results over the six catchments. Lastly, in Section 4, we summarize the results and present conclusions.

Study Area
Australia was selected as the study area to develop a simple scaling approach for LAI parameterization in hydrologic models. The region (10 • S-45 • S, 110 • E-155 • E) covers various climate and land cover types allowing to assess the general applicability of the proposed method, and its potential applicability elsewhere. As shown in Figure 1a,b, Australia consists of six primary land covers including forest, shrublands, woodlands, grasslands, croplands and unvegetated regions, and three dominant climate zones; tropical, arid and temperate. The presented land covers and climate zones are based on the moderate resolution imaging spectroradiometer (MODIS) yearly land cover data [32], and the updated Köppen-Geiger climate classification [33], respectively (see Section 2.2 for details).
For the simulations using SMART, as presented in Figure 1, we selected six ecohydrologically contrasted catchments from the hydrologic reference stations (HRS). These catchments are distributed over the mountainous regions and present non-stationary hydrologic responses [7]. The HRS consists of 222 monitoring sites with long-term streamflow records, and minimal flow regulation and land use change [34,35]. The catchment areas of HRS range from 101 to 232,846 km 2 , and the elevations vary from 284 to 1351 m (above sea level), respectively. Ajami et al. [7] identified that 20 out of 166 HRS catchments exhibit non-stationary hydrologic responses based on the long-term trends in annual runoff ratio, and subsequently developed an ecohydrologic catchment classification method using rainfall, runoff, vegetation productivity and evapotranspiration data to characterize ecohydrologic response. According to this classification, a positive or negative correlation between annual precipitation and fractional vegetation cover (productivity) in a catchment determines whether catchment-scale vegetation productivity is water limited (class A, positive correlation) or nutrient/light limited (class B, negative correlation). Both A and B classes are further subdivided into two groups based on the relationship between mean annual ET and vegetation productivity. In A1 catchments, productivity depends on the dominance of structural control (increases in LAI), and annual ET increases with annual productivity increase. In A2 catchments, physiological control (decreases in stomatal conductance) is dominant, and annual ET decreases with increases in productivity. In group B catchments, productivity is likely limited by biogeochemical factors. Productivity of B1 catchments is limited by nutrients, and in B2 catchments light availability is a limiting factor. Note that none of the 20 non-stationary catchments belongs to the A2 class. Here, we selected six catchments out of the 20 non-stationary catchments for hydrologic modeling. The selection included two catchments from each class with different proportion of land cover types. Details of the six catchments are summarized in Table 1. All of these catchments represent upland/hilly terrains, creating a mix of overland and lateral subsurface flow. Max Ele = maximum elevation above sea level; P = annual mean precipitation; and PET = annual mean potential evapotranspiration. Note that the land cover proportions (%) are based on the MODIS land cover product (F = forest; W=woodlands; G = grasslands; C = croplands). Ajami et al. [7] identified that 20 out of 166 HRS catchments exhibit non-stationary hydrologic responses based on the long-term trends in annual runoff ratio, and subsequently developed an ecohydrologic catchment classification method using rainfall, runoff, vegetation productivity and evapotranspiration data to characterize ecohydrologic response. According to this classification, a positive or negative correlation between annual precipitation and fractional vegetation cover (productivity) in a catchment determines whether catchment-scale vegetation productivity is water limited (class A, positive correlation) or nutrient/light limited (class B, negative correlation). Both A and B classes are further subdivided into two groups based on the relationship between mean annual ET and vegetation productivity. In A1 catchments, productivity depends on the dominance of structural control (increases in LAI), and annual ET increases with annual productivity increase. In A2 catchments, physiological control (decreases in stomatal conductance) is dominant, and annual ET decreases with increases in productivity. In group B catchments, productivity is likely limited by biogeochemical factors. Productivity of B1 catchments is limited by nutrients, and in B2 catchments light availability is a limiting factor. Note that none of the 20 non-stationary catchments belongs to the A2 class. Here, we selected six catchments out of the 20 non-stationary catchments for hydrologic modeling. The selection included two catchments from each class with different proportion of land cover types. Details of the six catchments are summarized in Table 1. All of these catchments represent upland/hilly terrains, creating a mix of overland and lateral subsurface flow. Max Ele = maximum elevation above sea level; P = annual mean precipitation; and PET = annual mean potential evapotranspiration. Note that the land cover proportions (%) are based on the MODIS land cover product (F = forest; W=woodlands; G = grasslands; C = croplands).

Remote Sensing and Forcing Datasets
The datasets used in this study span over a 12-year study period from November 2002 to December 2014. As summarized in Table 2, these data include remotely sensed observations and gridded climate products and are grouped to three categories by their primary purposes: LAI modelling; Remote Sens. 2020, 12, 3051 5 of 19 SMART modelling; evaluation of SM, ET, and Q simulations. Although ground measurements are generally regarded as truth and used for model validation, availability of in situ data is considerably limited in time and space [36][37][38]. Therefore, we used satellite-derived data for the catchment-scale evaluations of SM and ET similar to previous investigations [39,40]. For the LAI modelling, the 12-year period was separated into two segments: 2002-2009 (seven years) for model calibration, and 2010-2014 (five years) for model validation, respectively. The SMART modelling was implemented over a period from 1 January 2010 to 31 December 2014, and the first-year of model simulation was excluded from the evaluation in order to minimize the impact of model initialization.
In this study, three remote sensing datasets were used for LAI parameterization and hydrologic model evaluation at catchment scales. MODIS-derived 8-day composite LAI (MCD15A2H from Terra and Aqua satellites, V006) products [41] were used to cover the 12-year study period. As the products consist of 10 • × 10 • tiles in the sinusoidal projection, they were re-projected and resampled to 0.05 • × 0.05 • grids (identical to the gridded rainfall data) by using the MODIS re-projection tool [53]. The Savitzky-Golay filter [54] in the TIMESAT software package [55] was used for removing outliers from the reprojected LAI time series [56][57][58]. The filtered 8-day LAI maps were linearly interpolated to a daily time scale for temporally matching with other hydro-climatological datasets (i.e., rainfall, potential/actual evapotranspiration, and surface soil moisture) used in this study.
The European Space Agency Climate Change Initiative (ESA CCI) surface SM data was used for validating the SMART-simulated SM at the topsoil layer (<10 cm). The ESA CCI has released daily surface SM (<2 cm) products gridded at 0.25 • that are derived from seven passive and three active microwave spaceborne instruments covering 39 years from 1 November, 1978 to 30 June, 2018 (V04.4) [49][50][51]. Among the three surface SM products: active, passive, and active-passive combined data, the active-passive combined data were used for model validation in this study. MODIS-derived 8-day evapotranspiration data (MYD16A2, from Aqua satellite, V600) were used for evaluating SMART simulated ET. The raw ET data consisted of 10 • × 10 • tiles in the sinusoidal projection and Remote Sens. 2020, 12, 3051 6 of 19 were resampled to 0.05 • × 0.05 • grids through the MODIS re-projection tool. Only good quality data were included by using MODIS pixel-level quality code data. The spatiotemporal gaps in the quality-assured ET data were simultaneously filled by using a generic smooth parameter (s = 10 −6 ) in the three-dimensional (3-D) gap-filling method, called penalized least-square regression based on 3-D discrete cosine transform (DCT-PLS) [59][60][61].
For the LAI and SMART modeling, the Australian water availability project (AWAP) rainfall data gridded at 0.05 • [42] were used. The AWAP product uses a thin plate smoothing spline to interpolate monthly rainfall climatology derived from ground measurements, and a successive correction method is used for interpolating anomalies of daily rainfall [62]. The AWAP rainfall data is available from 1900 onwards as operational daily maps over Australia with a day of latency, and has been frequently used for various hydro-climatological studies [63][64][65][66]. For hydrologic modeling, daily PET data gridded at 0.05 • from the Australian Water Resource Assessment Landscape (AWRA-L) were used [48]. This product uses the Penman equation [67] for PET estimation. Alternate satellite precipitation products covering global extents are becoming increasingly refined and available for future extensions of the study to other parts of the world [68].

Description of SMART
The Soil Moisture and Runoff simulation Toolkit (SMART) [28] was used for daily simulations over the six catchments described in in Section 2.1. SMART is a geographic information system (GIS)-based semi-distributed hydrologic modelling framework designed for large catchment-scale simulations. Computational efficiency in SMART is achieved by delineating topologically connected hydrologic response units (HRUs), and formulating a series of equivalent cross sections (ECSs) by weighting topographic and physiographic properties of a part or an entire first-order sub-basin. In the 3-ECSs delineation approach employed here [29,30], 2 or 3 ECSs are delineated for every sub-basin depending on the number of hillslopes in each sub-basin. Each ECS consists of 4 HRUs (upslope, mid-slope, foot-slope, and alluvial flats) delineated based on changes in slope versus distance from the river. In each HRU, the dominant landcover type is selected for hydrologic simulation. Previous application of the ECSs approach reported 3.7 to 22.8 times reduction in the number of computational units [30]. The current version of SMART adopts the Unsaturated Soil Moisture Movement Model (U3M-2D) to solve the 2-dimensional Richards equation for every ECS by using spatially distributed climate, land cover and soil type information. U3M-2D developed by Tuteja et al. [69] is an extended two-dimensional version of U3M-1D [31]. U3M-2D calculates the vertical water balance by solving the 1D Richard's equation and then uses the unsaturated form of Darcy's law to transfer horizontal flow down the hillslope. The model generates daily time series of horizontal flow, deep drainage, total runoff, overstory and understory transpiration, soil evaporation, and soil moisture at multiple depths for every HRU. The U3M-2D uses monthly LAI climatology and scales it based on monthly rainfall to partition PET into overstory and understory potential transpirations and soil evaporation. Scaled monthly LAI is calculated as: where LAI i, j and P i, j are the scaled monthly mean LAI and monthly mean rainfall, respectively, in ith month (i = 1, . . . , 12) of jth year; LAI i and P i are monthly mean LAI and precipitation climatology, respectively, in ith month (i = 1, . . . , 12), β is a scaling parameter with a default value of 1 in SMART [70]. Based on the LAI i, j of every land cover type, U3M-2D partitions potential evapotranspiration (PET) into overstory and understory potential transpirations (ET o and ET u ), and potential soil evaporation (ET g ) for three different land cover types (i.e., tree, crop, and pasture) as where L is a canopy light extinction coefficient (default value of 0.8 for trees and 0.6 for crop and pasture are adopted in SMART) [70].
The following datasets were used for the SM and ET simulations with SMART: (1) Shuttle Radar Topography Mission (SRTM)-derived 1-s (~30m) Hydrologically Enforced Digital Elevation Model (DEM-H) Version 1.0 [43] for representing the topographic and physiographic properties of the catchments and delineating HRUs; (2) 2012 MODIS-derived yearly land cover (LC) map gridded at 0.05 • (MCD12C1, V051) [32] to represent landcover types. As summarized in Table 3, the MCD12C1 product consists of the 17-class International Geosphere Biosphere Program (IGBP) classified to 6 primary classes: forest, shrublands, woodlands, grasslands, croplands, and unvegetated regions. Note that the primary land cover classes in the MODIS LC product were aggregated to the three LC types used in Equation (2) of SMART; (3) The updated Köppen-Geiger climate classification [71] to define climate zones in SMART, and apply different rainfall and potential evapotranspiration datasets for each zone; (4) Soil attribute maps of silt, clay, and sand proportions at multiple depths, and soil depth at a 3-s (≈90m) spatial resolution [44][45][46][47]; and (5) daily precipitation and PET data as described in Section 2.2.

Proposed LAI Parameterization Method
Tesemma et al. [25] tested various types of linear and nonlinear models to estimate monthly LAI by using multiple climate predictors over the Goulburn-Broken catchment in southeastern Australia. Their results showed that a sigmoid model based on rainfall and reference potential evapotranspiration performs best over the three different land covers: crop, pasture, and trees in the catchment. The sigmoid function has the form Remote Sens. 2020, 12, 3051 8 of 19 whereL AI is the estimated monthly LAI; α 1, 2, 3 are model parameters that should be calibrated for each land cover classes (i.e., crop, pasture and tree) by minimizing the root-mean-square error between the modelled and remotely sensed LAI over the region; P and PET are n-month moving average rainfall and reference crop evapotranspiration, respectively (n = 6 for crop and pasture, and n = 9 for tree).
Based on the results of Tesemma et al. [25], we propose a simplified one-parameter sigmoid scaling model to estimate daily LAI based on monthly mean LAI climatology and rainfall aŝ whereL AI t is the estimated daily LAI at the time of t; LAI i and P i are identical to those in Equation (1) and, in this study, were derived from the 12-year time series of MODIS LAI and AWAP rainfall respectively; k (≥0) is the scaling parameter calibrated based on remotely sensed LAI; P 30 is the mean daily precipitation during recent 30 days (i.e., from t-29 to t). The k parameter in Equation (4) represents sensitivity of LAI dynamics to rainfall at every grid cell, and its value was optimized by maximizing the Pearson correlation coefficient (R) between the MODIS LAI time series and predictedL AI t at each grid cell over the calibration period. The benefit of this LAI scaling is that LAI data can be estimated for regions or periods where high quality remote sensing data are not available.

Hydrologic Modeling Scenarios and Evaluation Metrics
Three approaches for incorporating LAI dynamics into SMART were examined in this study to assess their impacts on simulated ET, SM and Q. These scenarios are (1) Opt def : using the scaled monthly LAI of Equation (1) (default in SMART), (2) Opt MOD : directly using the MODIS LAI (see Section 2.2), and (3) Opt est : using the estimated daily LAI by Equation (4).
Three metrics were used for model evaluation. Those are Pearson correlation coefficients (R), Nash-Sutcliffe model efficiency coefficient (NSE) [72] and root mean squared error (RMSE) between the simulated and observed SM, ET, and Q. Considering systematic differences in the amplitude and mean values of satellite-derived datasets resulting from limited parameterization and imperfect retrieval algorithms [36], R would be a more reliable metric because it is less sensitive to the bias or amplitude of variations [73].
Note that both satellite-derived and estimated SM and ET data within a catchment boundary were spatially averaged on a daily time step for catchment-scale comparisons given the coarse resolution of satellite data products and sizes of the study catchments. It should be also noted that the current version of SMART does not have a routing module. Therefore, total runoff (horizontal flow + deep drainage) values from each ECS were simply area-weighted averaged to represent Q at the outlet, and this limitation should be considered in interpreting the results of simulated Q.

Results of LAI Modelling
We compared the two LAI estimation approaches across Australia using the Opt est and Opt def methods for the five-year validation period from 2010 to 2014. For the LAI estimation by Opt est , the parameter k of Equation (4) was calibrated over the seven-year calibration period from 2002 to 2009 (Figure 2a). Then LAI time series at each grid cell were estimated by Equation (4) using the calibrated k. Likewise, for the case with Opt def , time series of LAI for every grid were generated by Equation (1). For both approaches the correlation coefficient between the estimated and MODIS LAI data at each grid cell was calculated to assess accuracy of the LAI parameterization approach.  shows boxplots comparing the two sets of correlation coefficients (R) derived from Opt def and Opt def approaches conditioned by the six primary land cover classes described in Table 3.
methods for the five-year validation period from 2010 to 2014. For the LAI estimation by Opt est , the parameter k of Equation (4) was calibrated over the seven-year calibration period from 2002 to 2009 (Figure 2a). Then LAI time series at each grid cell were estimated by Equation (4) using the calibrated k. Likewise, for the case with Opt def , time series of LAI for every grid were generated by Equation Error! Reference source not found.. For both approaches the correlation coefficient between the estimated and MODIS LAI data at each grid cell was calculated to assess accuracy of the LAI parameterization approach. Figure 2b shows boxplots comparing the two sets of correlation coefficients (R) derived from Opt def and Opt def approaches conditioned by the six primary land cover classes described in Table 3. As shown in Figure 2b, R values have been considerably improved for Opt est compared to Opt def indicating better performance of the sigmoid function in capturing LAI dynamics compared to the default LAI scaling option of Equation (1). When the results are grouped by the MODIS-derived six primary land covers (Figure 2b), the LAI estimation method best performed in woodlands, grasslands, and croplands compared to forest, shrublands, and unvegetated regions (Table 4). Larger k values coincide with the arid regions indicating larger sensitivity of vegetation cover to rainfall amounts over semi-arid ecosystems covered by grasslands and croplands [74][75][76][77]. As presented in Table 4, relatively small k values were observed for forests, shrublands, woodlands and unvegetated regions, meaning that these vegetation types are less sensitive to the rainfall amounts.
At the catchment scale, Opt est resulted in considerable improvements in R values of simulated LAI and remotely sensed LAI compared to Opt def as summarized in Table 5. In general, estimated catchment average k values in study catchments are small and less than 0.10 as the dominant land cover is forest and woodlands. In catchment 405238, k is 0.13 because this catchment is mainly covered by grass and croplands.

Results of Catchment-Scale Simulations of SM, ET and Q
Note that Results of the modeling with SMART for the six catchments are summarized in Table 6. First, Opt MOD improves SM and ET simulations compared to those from Opt def in terms of all metrics (i.e., R, NSE, and RMSE), with results being almost identical to that of Opt est . Therefore, the proposed LAI estimation method has considerable potential for adoption in hydrological modelling. Comparisons of SM simulations across multiple catchments and LAI options indicate that, A1 (water limited) catchments tend to have the highest R values (>0.8) between simulated and observed SM followed by the B2 (light limited) and B1 (nutrient limited) catchments. However, no difference among various LAI options were observed in simulated SM of A1 catchments, and the Opt est approach improved simulated SM in light/nutrient limited catchments. In terms of simulated ET, relatively low R values are also observed in B1 compared to B2 and A1 catchments. While the improvement in R for Opt est against Opt def is 10% or more for B1 and B2 catchments, the two water-limited catchments (A1) do not show such improvements in Opt est against Opt def . One possible explanation for getting similar results in the A1 catchments is that ET in A1 catchments is mainly limited by the precipitation amount, i.e., water supply [7]. As the current version of SMART does not have a routing module, it would be useful to focus on the differences in the metrics across the vegetation parameterization methods particularly in large catchments. For example, simulated discharge in the largest catchment (405238) showed notably poor performance with R of around 0.2 and NSE of −0.54 (Table 6). Even considering such limits, the differences in the Q simulations metrics are not considerable among different LAI parameterizations as shown in the last three columns of Table 6. This result indicates that simulated Q in SMART is less sensitive to the vegetation dynamics compared to those of the SM and ET. To further investigate the results, time series of satellite-derived and simulated SM, ET and Q are presented in Figures 3-5. Note that for better interpretation and visualization, only 1-year (i.e., 2014) of simulated results were presented for the Opt est and Opt def scenarios. As shown in Figure 3a-f, Opt def parameterization frequently underestimates top layer soil moisture compared to Opt est , especially during the dry period from January to May, while similar performance during very wet periods are also observed. Underestimation of top layer soil moisture further decreases deep drainage from root zone in models parameterized with Opt def . For the ET simulation results (Figure 4a-f), underestimations of SM during dry periods for Opt def results in ET overestimation. Again, both Opt est and Opt def do not show significant differences in the water-limited catchments (A1) despite differences in vegetation types. However, the differences in performance are larger for the nutrient-limited catchments (B1). For the Q simulation results (Figure 5a-f), it should be noted that the area weighted total runoff from all sub-basins instead of routed runoff does not properly represent the observed Q with attenuated and transformed peaks. In particular, these distortions are most noticeable in flows from the largest catchment (410734), and the two water-limited catchments (405238 and 410061) where small to medium peaks during the recession periods are disappearing. Despite these limitations, differences in the Q simulations by Opt est and Opt def are detectable in the rise and recession stages of Q. Small to moderated flow events from Opt est tends to be higher than that of Opt def , which can be explained by the under-and overestimations of SM and ET by Opt def , respectively.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 22 that of Opt def , which can be explained by the under-and overestimations of SM and ET by Opt def , respectively. While general conclusions cannot be drawn regarding the impacts of LAI dynamics on simulated hydrologic fluxes in ecohydrologically contrasted catchments given the limited number of study catchments, our results indicate improvements in overall model performance across the board. Further investigation is needed to show whether these conclusions hold across other ecohydrologically contrasted catchments. It should be noted that the impacts of LAI dynamics on simulated model states and fluxes depend on model structure as well, which is not explored here.        While general conclusions cannot be drawn regarding the impacts of LAI dynamics on simulated hydrologic fluxes in ecohydrologically contrasted catchments given the limited number of study catchments, our results indicate improvements in overall model performance across the board. Further investigation is needed to show whether these conclusions hold across other ecohydrologically contrasted catchments. It should be noted that the impacts of LAI dynamics on simulated model states and fluxes depend on model structure as well, which is not explored here.

Saptial Variability of SM and ET
In this section, we further compared the spatial variability of SM and ET derived from Opt est and Opt def at 318076 catchment. The reasons for selecting 318076 catchment are that the catchment size is relatively large (379.8 km 2 ) and improvements in simulated SM and ET are considerable compared to other catchments. As presented in Figure 6a, the elevation of 318076 catchment varies from 300 to 1500 m. Most of the land cover of this catchment consists of forest (88%) except for the southeastern highland area covered by grasslands (Table 1), and has considerably lower mean LAI than other regions (Figure 6b). Because of this spatial heterogeneity, different SM and ET values over the region are generally expected.

Saptial Variability of SM and ET
In this section, we further compared the spatial variability of SM and ET derived from Opt est and Opt def at 318076 catchment. The reasons for selecting 318076 catchment are that the catchment size is relatively large (379.8 km 2 ) and improvements in simulated SM and ET are considerable compared to other catchments. As presented in Figure 6a, the elevation of 318076 catchment varies from 300 to 1500 m. Most of the land cover of this catchment consists of forest (88%) except for the southeastern highland area covered by grasslands (Table 1), and has considerably lower mean LAI than other regions (Figure 6b). Because of this spatial heterogeneity, different SM and ET values over the region are generally expected.

Saptial Variability of SM and ET
In this section, we further compared the spatial variability of SM and ET derived from Opt est and Opt def at 318076 catchment. The reasons for selecting 318076 catchment are that the catchment size is relatively large (379.8 km 2 ) and improvements in simulated SM and ET are considerable compared to other catchments. As presented in Figure 6a, the elevation of 318076 catchment varies from 300 to 1500 m. Most of the land cover of this catchment consists of forest (88%) except for the southeastern highland area covered by grasslands (Table 1), and has considerably lower mean LAI than other regions (Figure 6b). Because of this spatial heterogeneity, different SM and ET values over the region are generally expected.   First, we performed two-sided paired t-tests at a significance level of 0.05 for the four pairs in Figure 7., i.e., a-b, c-d, e-f, and g-h, with the null hypothesis of the same mean for each pair. Results indicate that the mean values are significantly different (p-value < 0.05) for all pairs. In terms of the spatial patterns of SM and ET simulated by the two options, the differences in the dry day are larger than the wet day. This result is well explained by the spatial correlation coefficient of SM (ET) on the dry day which is 0.16 (0.77) and in the wet day is 0.89 (0.96). There are considerable differences in the mean and spatial pattern of SM during the dry day (Figure 7a,b) and the SM underestimation is explained in Section 3.2 and Figure 3.
This underestimated SM by Opt def can be explained by low LAI values estimated by Equation (1) over the two-month period from February to March (shaded area in Figure 8) while the MODIS LAI smoothly changes over that period. Lower LAI values in the default vegetation parameterization method (Opt def ) results in drier soil layer as shown again in November. The proposed method (Opt est ) compensates for these shortcomings and produces improved simulation results. First, we performed two-sided paired t-tests at a significance level of 0.05 for the four pairs in Figure 7., i.e., a-b, c-d, e-f, and g-h, with the null hypothesis of the same mean for each pair. Results indicate that the mean values are significantly different (p-value < 0.05) for all pairs. In terms of the spatial patterns of SM and ET simulated by the two options, the differences in the dry day are larger than the wet day. This result is well explained by the spatial correlation coefficient of SM (ET) on the dry day which is 0.16 (0.77) and in the wet day is 0.89 (0.96). There are considerable differences in the mean and spatial pattern of SM during the dry day (Figure 7a,b) and the SM underestimation is explained in Section 3.2 and Figure 3.
This underestimated SM by Opt def can be explained by low LAI values estimated by Equation Error! Reference source not found. over the two-month period from February to March (shaded area in Figure 8) while the MODIS LAI smoothly changes over that period. Lower LAI values in the default vegetation parameterization method (Opt def ) results in drier soil layer as shown again in November. The proposed method ( Opt est ) compensates for these shortcomings and produces improved simulation results.

Caveats and Follow-up Studies
In spite of the simplicity in the LAI modeling and improvements in the catchment-scale simulations, it is worth addressing the following improvements for practical applications of this framework in any hydrologic model.
First, the global monthly mean leaf area index climatology derived from more than three decades (1981-2015) of bi-weekly LAI data from the advanced very-high -resolution radiometer (AVHRR) gridded at 0.25° [78] could be an alternative dataset to the 12-year MODIS LAI data used herein for parameter estimation and obtaining the monthly mean LAI climatology. Second, we presented optimized k values for the six primary MODIS land covers as shown in Figure 2e. However, large uncertainties exist for k values for shrub, grass, and croplands indicating consideration of a more detailed land cover map may reduce uncertainty of estimated parameters for hydrologic modeling. Third, although we have shown simulated discharge from the different vegetation parameterization approaches, the results from large catchments raise the need for developing a streamflow routing module in SMART. Fourth, given the uncertainty of the reference datasets, it is necessary to use multiple data sources for a more comprehensive validation of ET and SM simulations [29,79,80]. Lastly, further improvements can be achieved by increasing the complexity of the LAI estimation model to include more climate variables, calibration parameters and/or use other model structures.

Caveats and Follow-Up Studies
In spite of the simplicity in the LAI modeling and improvements in the catchment-scale simulations, it is worth addressing the following improvements for practical applications of this framework in any hydrologic model.
First, the global monthly mean leaf area index climatology derived from more than three decades (1981-2015) of bi-weekly LAI data from the advanced very-high -resolution radiometer (AVHRR) gridded at 0.25 • [78] could be an alternative dataset to the 12-year MODIS LAI data used herein for parameter estimation and obtaining the monthly mean LAI climatology. Second, we presented optimized k values for the six primary MODIS land covers as shown in Figure 2e. However, large uncertainties exist for k values for shrub, grass, and croplands indicating consideration of a more detailed land cover map may reduce uncertainty of estimated parameters for hydrologic modeling. Third, although we have shown simulated discharge from the different vegetation parameterization approaches, the results from large catchments raise the need for developing a streamflow routing module in SMART. Fourth, given the uncertainty of the reference datasets, it is necessary to use multiple data sources for a more comprehensive validation of ET and SM simulations [29,79,80]. Lastly, further improvements can be achieved by increasing the complexity of the LAI estimation model to include more climate variables, calibration parameters and/or use other model structures.

Conclusions
In this study, a method was proposed for effectively incorporating LAI dynamics in a hydrologic model, SMART, by only using monthly mean LAI climatology and rainfall data over a region. The presented method adopts a simple sigmoid model that scales the monthly mean LAI climatology to a lower/higher value using a temporally aggregated rainfall amount and a parameter called k that control LAI-rainfall sensitivity. Here, the k parameter was optimized by maximizing the Pearson correlation coefficient between the estimated and satellite-derived LAI time series. It was found that forest, woodlands, and unvegetated regions have smaller k values, and present less sensitivity to the rainfall amounts compared to shrublands, grasslands and croplands that are more sensitive.
Surface soil moisture and evapotranspiration simulations illustrated that the proposed LAI estimation method performs almost identical to the case where the satellite-derived LAI time series over the six eco-hydrologically contrasted catchments in Australia are used. In nutrient/light limited catchments (B1 and B2), generally, correlation coefficients between simulated and observed SM and ET are improved by more than 10% compared to the default option in SMART. However, the choice of the LAI model did not impact SM and ET simulations in A1 catchments where vegetation productivity is water-limited. The differences in the discharge simulations by the default and proposed methods are noticeable during the rise and recession stages of hydrographs. Small to moderated flowrates from the proposed method are larger than that of the default method, which can be explained by the under-and overestimations of SM and ET by the default method during dry periods. The proposed LAI parametrization method compensates for the shortcomings of the default method which underestimates LAI and SM during dry periods. This can improve simulations especially in future climates given documented shifts towards more arid landscapes and drier pre-storm antecedent conditions that have been noted in the literature [81,82].
The LAI parameter estimation and model evaluations results showed promise, with possibilities of extension to other hydrologic models that need to include vegetation dynamics for periods or regions where satellite data are unavailable. Funding: This work has been undertaken as part of a Linkage Project (LP160100620) "Adapting catchment monitoring and potable water treatment to climate change" and we appreciate the funders, Australian Research Council, WaterNSW and Sydney Water.