Development and evaluation of L-THIA sub-daily flow and water quality prediction model Open Access

Hydrological changes due to urbanization and climate change are critically important in land use planning and water resources management, necessitating their quantification (Semadeni-Davies et al. 2008; Oudin et al. 2018). The influence of climate change on runoff can be directly estimated from annual climate and runoff time series based on the concept of runoff climate elasticity, or estimated through hydrological modeling (Sankarasubramanian et al. 2001; Chiew 2006; Chiew et al. 2009).

Recently, artificial intelligence (AI) techniques, such as machine learning and deep learning, have been widely employed to predict runoff and water quality (WQ) in watersheds worldwide (Parsaie & Haghiabi 2017; Qishlaqi et al. 2017; Haghiabi et al. 2018; Zhu et al. 2022; Yifru et al. 2024a, 2024b). However, hydrological processes are highly complex, involving the interaction of various factors such as rainfall, land use, soil properties, and human activities (McDonnell et al. 2007). These factors vary across time and space, making it challenging for data-driven models alone to fully capture these actions. While machine learning (ML) models have potential, they are limited by their black-box nature, practicality in watershed management, and difficulty in explaining physical hydrological processes (Yaseen 2023).

Physics-based watershed models, grounded in fundamental hydrological principles, offer the advantage of accurately simulating runoff and WQ changes under diverse environmental conditions. Such models can analyze various scenarios, including climate change, land use shifts, and the implementation of best management practices, providing essential support for long-term WQ management and decision-making.

The long-term hydrologic impact assessment (L-THIA) model was initially developed to estimate runoff, groundwater, and non-point source pollution resulting from land use changes (Bhaduri et al. 2001; Ma 2004). A key strength of L-THIA is its ability to estimate average annual runoff based on long-term climate data across different land use scenarios. However, the National Resources Conservation Service-curve number (CN) method in L-THIA has limitations in accurately calculating direct runoff during low-flow periods, as the static CN value fails to capture the relationship between rainfall and the asymptotic CN, particularly in cases of lower rainfall (Kwak et al. 2010; Ryu 2016). Thus, originally coded in Fortran, the model has been improved by several researchers to enhance runoff and WQ prediction accuracy. Ryu (2016) developed and evaluated the L-THIA ACN-WQ 2016 model, utilizing long-term flow-quality monitoring data observed in Korea to create 13 asymptotic CN regression formulas for different land cover types and an additional 52 formulas for various land cover and hydrologic soil group combinations, improving pollutant load evaluation. Kum (2018) further enhanced the L-THIA ACN-WQ 2016 model by developing an auto-calibration tool for optimizing flow and WQ parameters, accounting for hydraulic structures, and deriving a formula for estimating initial infiltration losses. However, the L-THIA ACN-WQ 2018 model remains limited to daily simulations due to its reliance on the soil conservation service (SCS)-CN method.

More than 70% of Korea's land is mountainous, with an average slope of 20%, causing rapid runoff of rainwater to rivers and a large coefficient of flow variation. Thus, sub-daily simulations are essential to capture both the dynamics of rainfall–runoff and the transport of associated pollutants (Jeong et al. 2010). Given the increasing variability in rainfall patterns due to climate change, models operating on shorter time scales tend to provide more reliable results (Ficchì et al. 2016). Studies have shown that sub-daily models outperform daily models in predicting direct runoff (Yu et al. 1998; Socolofsky et al. 2001; Kandel et al. 2005).

Maharjan et al. (2013) demonstrated that finer temporal resolution of rainfall data improves model performance, while Yang et al. (2016) found that sub-daily precipitation data yielded better peak flow predictions during the flood season than daily data in the upper Huai River watershed using the Soil and Water Assessment Tool (Arnold et al. 1998).

Therefore, a watershed model capable of long-term continuous and event-based simulations is needed to assess the impacts of urbanization and climate change. However, the current L-THIA model does not support hourly flow and WQ simulations at the watershed scale, nor does it integrate modules for direct and baseflow, channel routing, and WQ simulations in an open-source framework for sub-daily hydrological modeling. To address these limitations, this study aims to develop and evaluate a Python-based L-THIA sub-daily flow/WQ model for simulating hourly flow and WQ in agricultural and urban watersheds.

The previous L-THIA ACN-WQ models (Ryu 2016; Kum 2018) were developed in Fortran, which posed limitations, such as being difficult to modify. Therefore, this study developed the L-THIA sub-daily flow/WQ model based on Python. The model requires input data such as continuous rainfall data, hydrologic response unit (HRU) classified by sub-watershed, terrain data, land use map, and land cover map.

In a watershed with a high urbanization rate, runoff analysis should be performed using urban rainfall data, and independent rainfall events should be calculated using interevent time definition (IETD) (Lee & Chung 2017). Huff (1967) investigated the time distribution of 261 rainfall events over 12 years (1955–1966) using 49 rain gauge records on the 400 mi² network in eastern-central Illinois. In this investigation, rainfall was defined as an interval of 6 h or more from the previous and subsequent rainfall. Furthermore, the minimum preceding drying time suggested in the non-point pollution load evaluation technique study (National Institute of Environmental Research 2006) was also determined to be 6 h. In this study, rainfall events were separated by setting the IETD to 6 h, and a module for event separation was developed using continuous rainfall data.

This study has developed a direct runoff calculation module of the L-THIA sub-daily flow/WQ model using the CN4GA method suggested by Grimaldi et al. (2013). The CN4GA method is a combination of the Green & Ampt (1911) method as modified by Mein and Larson, known as the Green-Ampt Mein-Larson (GAML) model (Mein & Larson 1973) and SCS-CN methods. The total amount of net rainfall events and the initial loss calculated in the SCS-CN method are used to correct the ponding time and effective saturated hydraulic conductivity (K_eff) of the GAML method (Grimaldi et al. 2013). As many authors have emphasized, since this effective saturated hydraulic conductivity K_eff is very difficult to estimate, the value suggested in the literature is often used, so it is difficult to calculate the infiltration amount and direct runoff amount by reflecting the characteristics of the site.

The CN4GA method, developed by Grimaldi et al. (2013), improves surface runoff estimation at sub-daily time scales by combining the SCS-CN method with the Green–Ampt method. While the SCS-CN method calculates total runoff volume, it lacks temporal distribution for short-duration events. CN4GA addresses this by using the Green–Ampt model to distribute excess rainfall over time, allowing for more accurate runoff predictions, particularly in ungauged basins and extreme rainfall events.

Therefore, accurately estimating infiltration and the effective saturated hydraulic conductivity (K_eff) is critical when using CN4GA for hourly runoff calculations. In this study, the optimal K_eff in CN4GA is determined by adjusting the effective saturated hydraulic conductivity in the GAML model until the cumulative infiltration (F_GA) of the GAML model matches the cumulative infiltration (F_CN) calculated by the SCS-CN method. For this, FCN is calculated using the asymptotic CN method, which is based on land cover-specific monitoring data provided by the Ministry of Environment (MOE), applying the regression formulas derived by Ryu (2016).

Initial parameters for GAML, such as K_GA and cumulative infiltration F_GA, are derived from literature values of soil properties. If the calculated F_GA is greater than F_CN, K_eff is lowered to recalculate GAML's cumulative infiltration. Conversely, if F_GA is less than F_CN, K_eff is increased. This iterative process continues until the difference between F_GA and F_CN is less than 0.1, determining the optimal effective saturated hydraulic conductivity (K_eff-opt). This K_eff-opt value is then used to compute direct runoff and infiltration at each time step in the GAML model.

In this approach, F_CN is determined using an asymptotic CN based on the magnitude of the rainfall event, rather than a single CN calculated from land use and hydrological soil group combinations. Additionally, to ensure K_eff-opt reflects soil's physical properties, the minimum and maximum K_sat values from more than 19,000 soil databases, as proposed by García-Gutiérrez et al. (2018), were applied (Table 1). This procedure was repeated for all HRUs in the watershed and for each rainfall event calculated by the rainfall event calculation module.

Using HRU per filtration determined by the CN4GA method, the amount of water that flows into the aquifer and the amount of base runoff that flows into the river can be calculated using the the following equation:

Flood routing in rivers can be divided into hydrological flood routing, which interprets the storage equation, and hydraulic flood routing, which is solved by applying boundary conditions to partial differential continuity equations and motion equations of unsteady non-uniform flow. The hydrological flood routing method is simpler than the hydraulic flood routing method.

Hydrologic flood routing can be divided into reservoir routing, channel routing, and watershed routing. To estimate the runoff from the watershed outlet in the large watershed that is divided into several sub-watersheds, the channel routing method is used. Generally, the most widely used Muskingum method for channel routing was developed by McCarthy (1939), which establishes a relationship between storage, inflow, and outflow.

In the L-THIA sub-daily flow/WQ model, WQ is calculated from the direct and base flow of each HRU. Also, the model uses pollutant load per HRU that considers the event mean concentration (EMC) data investigated in the single land cover of non-point long-term runoff monitoring of the environmental basic survey project of the MOE. The calculated pollutant loads for each HRU are combined for each sub-watershed, then the WQ concentration is calculated using the runoff per sub-watershed, and this WQ concentration is used as the input data for QUAL2E (Ryu 2016; Kum 2018). QUAL2E is a WQ simulation model designed to predict changes in stream WQ under various environmental conditions (Brown & Barnwell 1987).

The QUAL2E model necessitates the inclusion of organic nitrogen (Org-N), nitrate–nitrogen (NO₃–N), nitrite–nitrogen (NO₂–N), and ammonia–nitrogen (NH₃–N) for nitrogen simulation, and organic phosphorus (Org-P), and inorganic phosphorus (Inorganic-P) for phosphorus simulation. However, it is difficult to use the EMC data of the non-point long-term outflow monitoring used in this study as input data for the QUAL2E model because it is presented only as T-N and T-P of nitrogen and phosphorus. Therefore, in the L-THIA ACN-WQ 2016 model and the L-THIA ACN-WQ 2018 model, the user determined the ratios of T-N and T-P, so that T-N could be divided into Org-N, NO₃–N, NO₂–N, and NH₃–N, and T-P could be divided into Org-P and Inorganic-P (Ryu 2016; Kum 2018).

Although the L-THIA ACN-WQ 2016 model and L-THIA ACN-WQ 2018 model have been developed to calculate various pollutant loads, there is a limitation that SS in the river cannot be calculated. Among non-point pollutants, soil and sand account for a large portion of rainfall runoff and cause various problems, such as adversely affecting WQ and aquatic ecosystems in downstream water systems.

Therefore, this study developed and integrated an SS calculation module into the L-THIA sub-daily flow/WQ model to evaluate suspended solids (SS) generated by rainfall runoff. To calculate the pollutant load per HRU, EMC data classified by land use and rainfall class from the MOE were used. SS transport in a channel network consists of two process functions, deposition, and degradation, which are calculated simultaneously and consist of a landscape component and a channel component.

In this study, the landscape component was calculated as the pollution load using the EMC data for each land cover, and the channel erosion component was calculated using the simplified Bagnold equation (Bagnold 1977), which estimates sediment transport based on flow velocity and stream power.

For the Bokha River, flow and WQ are measured every 8 days at Sangbaek Bridge, while sub-daily flow measurements are taken at Heungcheondaegyo Bridge. For the Gap River, flow and WQ are measured every 8 days at Bulmu Bridge, while hourly flow measurements are taken at Singu Bridge.

The Bokha River is a major tributary of the Namhan River, with non-point pollution sources such as farmlands scattered throughout the watershed (Nam et al. 2017). Additionally, discharges from various wastewater treatment plants in the Bokha River watershed have a significant impact on pollutant loads during the low-flow season (Lee et al. 2018b). Thus, the MOE manages the WQ of the Bokha River by setting target concentrations of biochemical oxygen demand (BOD₅) at 2.1 mg/L and T-P at 0.145 mg/L.

In the Gap River watershed, the WQ deteriorates due to the confluence of discharged water from the Daejeon sewage treatment plant located downstream, which is the cause of the deterioration of the WQ of the Geum River water system (Lee & Seo 2015). Additionally, the Gap River watershed experiences ongoing land use changes, such as urbanization, and shifts in precipitation patterns due to climate change, making it a priority area for the development of effective water pollution management strategies (Kwon et al. 2016; Yifru et al. 2024a). According to the National Institute of Environmental Research (2006), the pollutants targeted for management in the Geum River system from 2021 to 2030 include BOD₅ and T-P. For the Gap River A branch, the WQ targets are set at 4.1 mg/L for BOD₅ and 0.118 mg/L for T-P.

The input data of the L-THIA sub-daily flow/WQ model requires land use, soil, terrain data, and rainfall data. Among the input data of the model, the digital elevation model was created and used as topographical data using the digital topographic map provided by the National Geographic Information Institute (http://www.ngii.go.kr, accessed on 12 October 2024), for the land use map, the mid-class land use map in the 2010s provided by the Environment Spatial Information Service (https://egis.me.go.kr, accessed on 12 October 2024) of the MOE was used in this study. As for the soil map, the schematic soil map provided by the Rural Development Administration was collected and used. Baeksa and Daejeon meteorological station data were used for the hourly precipitation data in the Bokha River and Gap River watersheds, respectively, and were collected from the Korea Meteorological Administration.

To evaluate the hourly flow simulation performance of the L-THIA sub-daily flow/WQ model, the model was evaluated in the Bokha River and Gap River watersheds. Model parameters were calibrated using hourly flow data from the Bokha River (2020) and the Gap River (2018). For WQ, since hourly observations are not currently available in either watershed, calibration for flow, T-P, and SS was conducted using observed data measured every 8 days at Sangbaek Bridge for the Bokha River and Bulmu Bridge for the Gap River. During this process, discharge, T-P, and SS data from the wastewater treatment plants within each watershed were incorporated as boundary conditions (Table 2). Finally, the calibrated model was used to simulate hourly T-P and SS, and the results were aggregated into daily intervals to evaluate the model's predictive performance.

In this study, the coefficient of determination (R²), Nash–Sutcliffe efficiency (NSE), and percent bias (PBIAS) were used, as these are commonly applied evaluation metrics for such models. R² ranges from −1 to 1 and indicates the degree of linear relationship between observed and simulated data. NSE, which also ranges from −∞ to 1, measures the model's predictive accuracy, where a value closer to 1 indicates better performance. PBIAS describes the average tendency to be overestimated or underestimated by comparing simulated and observed values.

According to Moriasi et al. (2015), R² values between 0.6 and 0.75 are considered ‘Satisfactory’ values between 0.75 and 0.85 are rated as ‘Good’ and values above 0.85 are evaluated as ‘Very Good’. For NSE, values between 0.5 and 0.7 are classified as ‘Satisfactory’ values between 0.7 and 0.8 as ‘Good’ and values exceeding 0.8 as ‘Very Good’. PBIAS is evaluated based on its absolute value; a range of 10–15 is deemed ‘Satisfactory’ 5–10 is considered ‘Good’ and values below 5 are rated as ‘Very Good’.

A total of 52 rainfall events occurred at the Daejeon meteorological station, with a total of 1,533.2 mm of rainfall. The month with the most rainfall events was December, and a total of six rainfall events were calculated. This is considered to be due to the rainfall events with a duration of 42 and 48 h, in July and August. The rainfall event with the largest amount of rainfall and the longest duration at the Daejeon meteorological station was the rainfall event that lasted from 11:00 on 26 August 2018, to 11:00 on 28 August 2018. The highest rainfall intensity occurred on 28 August 2018, and lasted for 1 h at 6:00. In the case of the 29th rain event at the Daejeon meteorological station, when the no-rain time is applied to 6 h, it is calculated as one event, but when the no-rain time is reduced, it is judged to be calculated as 3–4 events.

When an outlier exists in the actual observed data during model calibration, it greatly affects the deduction and calibration of the optimal parameters. An outlier refers to data points that deviate substantially from typical observations (Jeon et al. 2020). Causes of outliers include data input errors, measurement errors, experimental errors, and sampling errors (Kim 2017). Various methods exist for detecting outliers, but in most time series data, the residual method is commonly applied. In the 2020 sub-daily flow data for the Bokha River watershed, there appears to be a section where runoff is excessively measured compared to rainfall. To identify these outliers, the standardized residual method is used, which detects outliers by examining the difference between the observed values and those predicted by the regression equation based on the sample data. When the standardized residual exceeds ±2, the data point is generally classified as an outlier.

The standardized residuals from 11:00 to 14:00 on 2 August 2020, were 16, 12, and −11, significantly exceeding the typical threshold for outliers, and were thus classified as outliers. Furthermore, an analysis of runoff rates for all events in the Bokha River watershed revealed that the average runoff rate for 2020 was 14%, with the runoff rate from July to September, during the summer, averaging around 20%.

According to Lee et al. (2018a), the average runoff rate of the Bokha River watershed was about 29%, which was the lowest among 27 points in the Han River water system. However, the runoff rate was about 76% from 11:00 on 2 August 2020, to 14:00 on 2 August 2020, when the sub-daily runoff was observed to be 1,403.1 m³/s. This is considered a measurement error, as the runoff amount was excessively high compared to the recorded rainfall. Based on this analysis, the period from 11:00 to 14:00 on 2 August 2020, was identified as an outlier. When this period was excluded, the runoff rate was found to be approximately 22%, which is not significantly different from the average summer runoff rate. Therefore, when calibrating the model, rainfall events with excessively observed runoff rates were classified as outliers and then removed, and river watershed calibration was conducted.

Therefore, in the calibration results of this study, both R² and NSE met the evaluation criteria. Notably, for PBIAS, the Bokha River watershed was rated as ‘Very Good’ and the Gap River watershed as ‘Good,’ indicating that the natural phenomenon was well simulated.

As a result, except for event 29, R² ranged from 0.66 to 0.96, and NSE ranged from 0.59 to 0.8. In the case of the 29th rainfall event, unlike the other events, rainfall exceeding 9 mm per hour occurred nine times, with intervals ranging from a minimum of 3 h to a maximum of 5 h. It is considered that multiple rainfall events appeared as one continuous event. Accordingly, the simulated flow in the latter part of the rainfall event appears to differ from the observed flow.

In the downstream region of the Bokha River watershed, rice farming is prevalent. Korean rice paddies are typically submerged during the irrigation season, which generally begins in spring, resulting in agricultural water entering the river and affecting both flow and WQ (Lee et al. 2023). Particularly, extensive groundwater use for water curtain cultivation increases groundwater extraction, leading to degraded river WQ during the dry season (Kang et al. 2022). However, the L-THIA sub-daily flow/WQ model did not account for the influence of groundwater and agricultural irrigation water, likely contributing to the underestimation of WQ during periods of low rainfall. The overestimation of WQ between November and January is likely due to the greater emphasis placed on high-concentration data from significant rainfall events during the summer, which were heavily reflected in the calibration process. Moreover, soil moisture is a key factor in controlling hydrological processes such as evaporation, infiltration, and runoff, making it crucial for accurate watershed hydrological modeling (Al-Yaari et al. 2019).

Thus, considering the physical behavior of soil moisture in the soil layer is essential for accurately predicting flow and WQ at sub-daily intervals (Qi et al. 2018). This requires precise estimation of spatially distributed soil parameters and sophisticated modeling of soil moisture (Shin et al. 2022; Lee et al. 2024). However, this L-THIA sub-daily flow/WQ model simulates soil moisture using a simplified storage routing method based on daily processes, which introduced parameter uncertainty, causing the model to overestimate WQ during low-flow periods when the influence of rainfall was minimal.

Overall, the calibration results demonstrate that the L-THIA sub-daily flow/WQ model performs reasonably well in simulating flow, WQ, and SS in both watersheds. However, the methods for estimating surface flow differ between hourly and daily hydrological simulations, and the sensitive parameters also vary (Jang & Kim 2016). To improve the estimation of more reasonable hydrological parameters for hourly simulations, securing WQ data for hourly model calibration is essential. Furthermore, additional research is needed to account for regional environmental factors such as groundwater influence, soil moisture, and agricultural practices to enhance model accuracy in future applications.

The purpose of this study was to develop and evaluate the L-THIA sub-daily flow/WQ model, which enables hourly flow simulations using the asymptotic CN regression equation and the CN4GA method. The model was applied to two distinct watersheds: the Bokha River watershed, with a high proportion of agricultural land, and the Gap River watershed, dominated by urban areas.

The results demonstrated that the L-THIA sub-daily flow/WQ model performs reasonably well in simulating both flow and WQ at sub-daily intervals in both watersheds. The model successfully captured peak runoff during rainfall events, achieving satisfactory performance for key flow and WQ indicators. However, some inaccuracies were observed during low-flow periods, likely due to the simplified method of soil moisture and the absence of groundwater and agricultural irrigation influences.

Through these results, we confirmed that the model's current limitations stem from the simplified soil moisture modeling and the exclusion of groundwater and irrigation water contributions, which are critical during periods of low rainfall. Improving the model to incorporate more sophisticated soil moisture dynamics and account for these additional factors will enhance its predictive accuracy, particularly for low-flow conditions. Moreover, the rainfall event calculation module, which currently uses a fixed 6-h minimum non-rainfall period, may need to be adjusted based on regional rainfall characteristics, requiring further research to optimize this parameter. Additionally, considering the diverse nonlinear watershed and climatic characteristics, utilizing advanced AI techniques for rainfall event detection could prove beneficial in improving accuracy and adaptability across different regions.

While the model has significant potential for application in diverse watersheds, there are uncertainties, especially in regions where complex interactions with groundwater exist, that should be carefully considered. Securing sufficient hourly observational data are also crucial for validating the model's performance and improving its reliability across various hydrological conditions. Additionally, the current model is limited by the lack of integration with robust calibration algorithms, such as SUFI-2 and genetic algorithms, making accurate parameter estimation challenging.

Despite these limitations, the L-THIA sub-daily flow/WQ model shows considerable potential for evaluating point and non-point source pollution in various watersheds. By utilizing asymptotic CN and EMC methods that reflect land cover-specific rainfall–runoff behavior, the model provides valuable insights into watershed discharge and pollutant loads, supporting the formulation of WQ management strategies. Additionally, the model was developed using Python, enhancing its scalability and ease of use, and making it accessible to non-experts when integrated with open-source platforms like QGIS. Overall, this model is expected to be effectively utilized for predicting runoff and assessing pollutant loads during rainfall events across different watersheds, contributing significantly to the development of watershed management strategies.

This research was supported by the Korea Environment Industry & Technology Institute (KEITI) through the Aquatic Ecosystem Conservation Research Program funded by the Korean Ministry of Environment (MOE), Grant Number 2020003030004 and the Ministry of Environment of Korea as The SS (Surface Soil conservation and management) projects [2019002820003].

All relevant data are included in the paper or its Supplementary Information.

The authors declare no conflict of interest.