Enhanced Watershed Modeling by Incorporating Remotely Sensed Evapotranspiration 1 and Leaf Area Index 2 3

Enhanced Watershed Modeling by Incorporating Remotely Sensed Evapotranspiration 1 and Leaf Area Index 2 3 Sangchul Lee a*, Dongho Kim a, Gregory W. McCarty b, Martha Anderson b, Feng Gao b, Fangni 4 Lei b, Glenn E. Moglen b, Xuesong Zhang b, Haw Yen c, Junyu Qi d, Wade Crow b, In-Young 5 Yeo e, Liang Sun f 6 7 a School of Environmental Engineering, University of Seoul, Dongdaemun-gu, Seoul 02504, 8 Republic of Korea 9 b USDA-ARS, Hydrology and Remote Sensing Laboratory, Beltsville, MD 20705, USA 10 c Crop Science, Bayer U.S., 700 W Chesterfield Pkwy W, Chesterfield, MO 63017, USA 11 d Earth System Science Interdisciplinary Center, University of Maryland, College Park, 5825 12 University Research Ct, College Park, MD 20740, USA 13 e School of Engineering, the University of Newcastle, Callaghan NSW 2308, Australia 14 f Key Laboratory of Agricultural Remote Sensing, Ministry of Agriculture / Institute of 14 15 Agricultural Resources and Regional Planning, Chinese Academy of Agricultural Sciences, 15 16 Beijing 100081, China 17


Introduction 48
One major concern with regard to any hydrological modeling exercise is predictive uncertainty. 49 Although the reliability of the simulated outcomes is assessed via model calibration and validation 50 to some degree, predictive uncertainty always exists (Arnold et  products in reducing the degree of equifinality, which is the tendency for different parameter sets 82 (referred to as PARs hereafter) to produce equally acceptable model outputs (Beven, 2006). A 83 study by Rajib et al. (2018) found substantial improvement in the modeled ET predictions by 84 including vegetation parameters and the utility of RS-ET products in evaluating ET variations 85 across a landscape, indicating a change in the model performance measure, that is the Kling-Gupta 86 Efficiency (KGE) from 0.6 to 0.7. Thus, access to RS-ET products enables the assessment of the 87 model capacity to predict the spatial distribution of hydrologic variables (Becker et  The primary goal of this study was to explore the usefulness of the two remotely sensed 105 datasets, namely RS-ET and RS-LAI, in enhancing watershed model uncertainty for a small 106 watershed (221 km 2 ) within the coastal plain of the Chesapeake Bay Watershed (CBW). The 107 hydrologic model chosen for this study was SWAT because remotely sensed data have been widely 108 incorporated into this model. To achieve this research goal, this study conducted a lumped 109 parameterization at the watershed level using three constraints: streamflow, RS-ET, and RS-LAI 110 products. The PARs that resulted in acceptable streamflow and ET simulations (referred to as 111 "PARs-1," hereafter) were taken from all PARs explored for calibration. In addition, the PARs 112 with acceptable model performance measures for streamflow, ET, and LAI (referred to as "PARs-113 2," hereafter) were extracted from all explored PARs. The specific objectives of this study were 114 to: (i) compare the two PARs (i.e., PARs-1 and PARs-2) along with their simulated outputs (e.g., 115 streamflow, ET, and LAI), and explore the role of vegetation constraints (i.e., RS-LAI products) 116 and the Atlantic Ocean, resulting in fairly uniform precipitation over the course of the year (Fisher 139 et al., 2010). The study site is characterized by flat topography (0 -32 m above sea level). Irrigation 140 for corn and soybean production during the summer season has seen a substantial increase in this 141 region (Wolman, 2008), which amplifies water loss by ET during the summer season. Water 142 balance cycling in this region is greatly affected by seasonal variations in ET. Thus, an accurate 143 ET simulation for this region is crucial for advancing predictions from hydrological models. water balance of each HRU is computed using Eq. 1: 164 (1) 165 Where, is the final soil water content (mm H2O), 0 is the initial soil water content (mm 166 Ampt infiltration method (Green and Ampt, 1911). The former was used in this study. 172 The SWAT model first calculates potential ET (PET) and then estimates actual ET (AET). This study used the Penman-Monteith method, which is expressed in Eq. (2), as follows: 182 Where, E is the latent heat of vaporization (MJ kg −1 ), the depth rate evaporation (mm d -1 ), Δ 184 the slope of the saturation vapor pressure-temperature curve (kPa •C −1 ), the net radiation 185 (MJ m −2 d −1 ), G the ground heat flux density (MJ m −2 d −1 ), the air density (kg m −3 ), the 186 specific heat at constant pressure (MJ kg −1 •C −1 ), 0 the saturation vapor pressure of air at height 187 z (kPa), e the water vapor pressure of air at height z (kPa), γ the psychrometric constant (kPa •C 188 −1 ), r c the plant canopy resistance (s m −1 ) and r a the diffusion resistance of the air layer 189 (aerodynamic resistance) (s m −1 ). 190 In SWAT, dynamic LAI estimates are generated as a function of the optimal leaf area 191 development curve. This curve controls LAI growth through accumulated potential heat units. A 192 daily potential heat unit is computed as the difference between the daily average temperature and 193 base temperature. The base temperature is the minimum temperature for vegetation growth, and 194 its default value is set to 0 °C. If the base temperature is greater than the daily average temperature, 195 the daily heat unit is zero. During the initial growth period, leaf area development is simulated as 196 a function of the optimal leaf area development curve. 197 Where, is the fraction of the plant's maximum leaf area index corresponding to a given 199 The SWAT model requires climate and geospatial data as inputs for simulations (Table 1) Table 1   Daily streamflow records from 2010 to 2014 were obtained from USGS gauging station 239 #01491500, located at the outlet of the TCW (Fig. 1a). Daily RS-ET products were generated from 240 the regional Atmosphere-Land Exchange Inverse (ALEXI) model (Anderson et al., 1997(Anderson et al., , 2007  The study watershed was divided into 19 sub-watersheds that ranged between 0.09 and 32 266 km 2 . In the HRU generation process, the threshold area values of land use, soil, and slope were set 267 to >10%, >15%, and >15%, respectively. There were 542 HRUs (312 cropland HRUs and 39 forest 268 because they accounted for more than 90% of the watershed. In addition, corn and soybean 292 parameters were adjusted because the distribution and rotation of the two crops were well captured 293 by the land use map used in this study. The detailed practice schedules (e.g., the application timing 294 and amount of fertilizer, planting, and harvesting timings) of the two crops were developed by 295 local experts (Lee et al., 2016). Thus, the growth dynamics of corn and soybean were depicted in 296 our simulations. The double crop soybean was not calibrated as all the information described above 297 was made for summer crops. Table 3 lists the calibrated parameters and allowable ranges. 298

Spatial evaluation at sub-watershed level 324
The simulated ET and LAI were compared with RS-ET and RS-LAI products at the sub-325 watershed level. The RS-ET and RS-LAI products were discretized by the sub-watershed boundary 326 generated from the ArcSWAT interface using the input DEM (Winchell et al., 2007). The TCW 327 included 19 sub-watersheds. Except for one sub-watershed that was smaller than the LAI pixel 328 size (0.25 km 2 ), 18 sub-watersheds were used for the sub-watershed-level spatial evaluation. This 329 evaluation was conducted using PARs-1 and PARs-2 simulations. Furthermore, the KGE values 330 were computed for ET and LAI for individual sub-watersheds and the median KGE values. The 331 PARs with median KGE values greater than 0.5 for both ET and LAI were considered to represent 332 acceptable performance measures for the spatial distribution of ET and LAI at the sub-watershed 333 level. PARs that did not meet these criteria were viewed as unable to capture the spatial distribution 334 of ET and LAI at the sub-watershed level, although they showed acceptable performance at the 335 watershed level. The evaluation results were used to further assess the degree of equifinality. 336 337

Impacts of vegetation data on ET predictions and predictive uncertainty at the 339 watershed level 340
The watershed-level calibration results show that there were 14 PARs-1 and 6 PARs-2 (Table  341   4 The degree of equifinality was reduced from 14 to 6 with the inclusion of the RS-LAI. 348 Although RS-LAI was incorporated, a 50% reduction in equifinality was observed because both 349 the ET calculation and RS-ET considered the LAI. The ET calculation method in this study 350 (Penman-Monteith) used canopy resistance as a key variable, which was calculated from the LAI 351 in SWAT (Neitsch et al., 2011). RS-LAI data were used as inputs for RS-ET retrievals (Sun et al., 352 2017 The observed streamflow, RS-ET, and RS-LAI were plotted against the simulation results from 360 PARs-2 (Fig. 2). The simulated streamflow did not capture the observed peak flows during the 361 simulation period. This may be because the precipitation data collected at the weather stations did 362 not fully represent the spatial variations in meteorological conditions across the entire study site. As compared to streamflow and RS-LAI, low KGE values were observed in the ET simulations 384 (Fig. 2). Low accuracy of ET in this study was likely attributable to the exclusion of irrigation 385 practices in our simulations because of inadequate associated information, whereas the thermal ET 386 remote sensing approach directly captured the impact of irrigation on ET (Hain et al., 2015). A 387 previous study found that improved ET simulation resulted from the inclusion of irrigation 388 practices in the simulations (Chen et al., 2017). Depressional wetlands, which are abundant in 389 forested areas of this region, are likely to lose water via ET at rates higher than those captured by 390 the SWAT model. Therefore, the ET module in the forested settings could have been an additional 391 factor that led to low KGE values of ET (Fig. 2). Simulated LAI values were mostly lower than 392 observations during the winter season (Fig. 2). Winter cover crops are widely implemented in this 393 region to reduce nutrient loads. These crops have been shown to increase the winter vegetation 394 index (Hively et al., 2020). The omission of winter cover crops from the simulation used in this 395 study resulted in a low simulated LAI during the summer season. 396 397

Comparing model results with RS-ET and RS-LAI at the sub-watershed level 398
Sub-watershed-level KGE values were calculated for daily ET and LAI, as shown in Fig. 3 With respect to the sub-watershed results, the number of acceptable PARs decreased from six 411 (PARs-2) to three, which suggested that the sub-watershed-level assessment helped identify the 412 PARs that satisfactorily characterized internal processes at a finer spatial level. This finding 413 supports the conclusion that spatial assessment using remotely sensed data can further narrow the 414 acceptable PARs, thus reducing predictive uncertainty (e.g., equifinality). our study site (220 km 2 ) was small. Moreover, this study focused on the use of multiple remotely 435 sensed datasets to reduce predictive uncertainty. Therefore, the lumped parameterization used in 436 this study was sufficient to assess the prediction accuracy of the spatial distributions of ET and 437 LAI. 438 439

Limitations and implications 440
This study aimed to improve model predictions by accommodating remotely sensed ET and 441 LAI in an effort to contribute to watershed modeling. However, this study had several limitations 442 to be conisdered for future studies. Remotely sensed data inevitably include uncertainties that are 443 greater than those in observations collected at the watershed outlet (Vervoort et al., 2014) but they 444 also enable hydrological models to be evaluated at a finer spatial level than watersheds (Rajib et 445 al., 2018). Thus, the uncertainty embedded in remotely sensed data must be carefully considered 446 when incorporating remotely sensed data into watershed modeling. Furthermore, simulated ET 447 and LAI are highly influenced by the climatic data. In this study, three sets of climatic input data 448 (i.e., humidity, solar radiation, and wind speed) were prepared using SWAT's built-in weather 449 generator. This has also been practiced in previous studies ( SWAT's built-in weather generator computes solar and relative humidity by a function of 512 precipitation and temperature. Solar radiation and relative humidity are determined based on the 513 number of dry or wet days per given month. Solar radiation is assumed to be lower on wet day 514 ( ) and that the wet day solar radiation is the half of the dry day solar radiation ( ). 515 (2) 517 Where, is the average daily solar radiation for the month, is the total number of days in 518 the month, and are the total number of wet and dry days in the month, respectively. 519 To incorporate the effect of clear and overcast weather on generated values of relative humidity, 520 monthly average relative humidity values can be adjusted for wet or dry conditions. The wet day 521 average relative humidity is assumed to be greater than the dry day relative humidity by some 522 fraction as Eq. (3). The dry day relative humidity is computed as shown in Eq. (4). (4) 526 Where, ℎ is the average relative humidity of the month on wet days, ℎ is the average 527 relative humidity of the month on dry days, is a scaling factor that controls the degree of 528 deviation in relative humidity caused by the presence or absence of precipitation, ℎ is the 529 average relative humidity for the month, and are the number of wet days in the 530 month and the total number of days in the month, respectively. 531 Wind speed is generated for the potential evapotranspiration by the Penman-Monteith equation. 532 Mean daily wind speed is generated using the equation below.