Evaluation of Different Objective Functions Used in the SUFI-2 Calibration Process of SWAT-CUP on Water Balance Analysis: A Case Study of the Pursat River Basin, Cambodia

Abstract

Many calibration techniques have been developed for the Soil and Water Assessment Tool (SWAT). Among them, the SWAT calibration and uncertainty program (SWAT-CUP) with sequential uncertainty fitting 2 (SUFI-2) algorithm is widely used and several objective functions have been implemented in its calibration process. In this study, eight different objective functions were used in a calibration of stream flow of the Pursat River Basin of Cambodia, a tropical monsoon and forested watershed, to examine their influences on the calibration results, parameter optimizations, and water resources estimations. As results, many objective functions performed better than satisfactory in calibrating the SWAT model. However, different objective functions defined different fitted values and sensitivity rank of the calibrated parameters, except Nash–Sutcliffe efficiency (NSE) and ratio of standard deviation of observations to root mean square error (RSR) which are equivalent and produced quite identical simulation results including parameter sensitivity and fitted parameter values, leading to the same water balance components and water yields estimations. As they generated reasonable fitted parameter values, either NSE or RSR gave better estimation results of annual average water yield and other water balance components such as annual average evapotranspiration, groundwater flow, surface runoff, and lateral flow according to the characteristics of the river basin and the results and data of previous studies. Moreover, either of them was also better in calibrating base flow, falling limb, and overall the entire flow phases of the hydrograph in this area.

1. Introduction

Several distributed hydrological models have been increasingly and widely employed, which are extremely essential to simulate hydrological phenomena, to manage and plan water resources, to investigate sedimentation and water quality, and to forecast the impact of universal climate and land-use changes. The Soil and Water Assessment Tool (SWAT) is one of the most well-known and has demonstrated its strengths in the above aspects with a large and growing number of model application in various studies ranging from catchment to continental scales [1,2,3,4,5]. It is open source software developed by the U.S. Department of Agriculture (USDA) Agricultural Research Service and Texas A&M. This physically based, watershed-scale, continuous model is commonly applied for measuring the impact of land-use and climate changes and for determining the different watershed management activities on hydrology, sedimentation, and water quality [6,7,8,9,10]. However, proper calibration/validation and uncertainty analysis are required to improve the model simulation and analysis for the watershed.

Calibration and uncertainty analysis are intimately linked, and no calibration results should be presented without quantification of the degree of uncertainty of the model prediction [11]. To run a physically based distributed hydrological model, such as SWAT, various parameters must be calibrated because of measurement difficulties. Calibration is performed by carefully selecting values with their respective uncertainty ranges and by comparing the simulation output for a given set of measured data [7]. Several calibration techniques have been developed for SWAT, including manual and automated calibration procedures using the shuffled complex evolution method and other common methods. The SWAT calibration and uncertainty program (SWAT-CUP) links sequential uncertainty fitting 2 (SUFI-2), generalized likelihood uncertainty estimation (GLUE), parameter solution (ParaSol), Markov chain Monte Carlo (MCMC), and particle swarm optimization (PSO) algorithms to SWAT [12] and has been widely used for calibration and uncertainty analysis [13,14,15].

Among these algorithms, SUFI-2 has gained the most popularity among users for carrying out parameterization, sensitivity analysis, calibration, validation, and uncertainty analysis of hydrological parameters [16]. This is due to the availability of many parameters regarding water balance modeling, the easier calibration process to carry out within realizable time bounds, better accountability for uncertainties, and the smallest number of computational runs to achieve good prediction uncertainty bands and model performance [17,18,19]. However, many objective functions have been incorporated into the SUFI-2 calibration procedure, but most studies have not considered the influence of the choice of objective functions on the calibration results, the calibrated parameter values/ranges, or the estimated water balance components. One study conducted in two Iranian watersheds by Kouchi et al. [20] was likely the only one thus far that has compared the choice of objective functions during the calibration process in SUFI-2. That study showed that all objective functions used in the calibration had a similar performance; however, each produced notably different parameter ranges and consequently led to significant differences in the water resource estimates.

Thus in this study, the objectives were to examine the influences of eight different objective functions in the calibration process of SWAT-CUP using the SUFI-2 algorithm on the calibration results, parameter optimizations, and water balance component estimations when simulating the stream flow of the tropical monsoon and forested watershed of the Pursat River Basin of Cambodia. Those objective functions include the coefficient of determination (R²), modified coefficient of determination (bR²), Nash–Sutcliffe efficiency (NSE), modified Nash–Sutcliffe efficiency (MNS), ratio of standard deviation of observations to root mean square error (RSR), ranked sum of squares (SSQR), Kling–Gupta efficiency (KGE), and percent bias (PBIAS). We also aimed to identify a proper objective function which produced reasonable calibration results, calibrated parameter sets, and water resources estimations which reflected the characteristics of this river basin.

2. Study Area

The Pursat River Basin, located in Pursat Province of Cambodia, is one of the Tonle Sap sub-basins and drains an area of 5955 km² [21] (Figure 1a). The river flows for approximately 150 km from the Southwest at the Cardamom Mountains to the Northeast direction until the Tonle Sap Lake with a basin elevation ranging between 6 and 1717 m above sea level. There is an operating hydrological station named Bak Trakuon with a drainage area of 4245 km² [22] where this study focused on (Figure 1b). More than 75% of the river basin encompasses a hilly terrain, with an elevation greater than 30 m above sea level and is covered by forested land of varying density, while the remaining low-lying land is occupied by agriculture [23]. Those over 75% of the forest cover are mainly concentrated in Bak Trakuon drainage boundary (Figure 1b) where the major soil types are Dystric Leptosols and Dystric Leptosols/Dystric Cambisols (LPd and LPd/CMd) and Gleyic Acrisols/Plinthic Acrisols (ACg/ACp) [22]. The climate in Pursat River Basin is influenced by tropic monsoon systems with distinct wet and dry seasons. The wet season, which extends from May to November, receives approximately 90% of the total annual rainfall. The dry season, which extends from December to April, is characterized by the prevalence of hot and dry air [22]. The rainfall in the area increases with elevation, but the annual totals vary considerably from year to year with the average ranging from 1200 to 1700 mm [24]. Daily maximum temperatures vary between 36 °C during the hottest months (April–May) and 32 °C during the coldest months (December–January). Daily minimum temperatures vary between 25 and 17 °C. The annual average temperature is approximately 28 °C. The monthly mean relative humidity ranges from 66% in the dry season to 71% in the wet season with a mean annual humidity of 70% [22].

3. Materials and Methods

Figure 2 shows the workflow diagram for this study. To simulate hydrological processes in the river basin using the SWAT model, input data, such as DEM, soil type, land-use map, and climate data, were required. The simulation output of the model without any calibration is called the default/initial result. The initial result was checked by using the SWAT Check program to identify any potential model problems, and the simulated discharge was compared with the observed data. If there was a problem or a huge difference between the simulated and observed data, the model inputs (e.g., spatial data, rainfall) need to be focused on, and if there were no problems, model calibration/uncertainty analysis was performed. SWAT-CUP with SUFI-2 algorithm was employed to conduct the model calibration/uncertainty analysis. This process started with parameterization in which several parameters were selected with indication of their value ranges. This was based on the difference between the simulated and observed data, watershed characteristics, and literature reviews. The next step involved applying different objective functions before beginning the calibration and validation. The results of the calibration and validation were evaluated. If the criteria was not satisfied, parameterization was rechecked; then, the model was recalibrated and revalidated. After the calibration and validation process using different objective functions, each objective function was evaluated based on its calibration result, calibrated parameters, the estimation of discharge processes, and the simulation of the hydrograph components. The details of each step are described in the following sections.

3.1. SWAT Input Datasets

Input data, such as topography, land-use map, soil type, and weather data, were collected from different sources and agencies, as listed in Table 1. The digital elevation model (DEM) was from the Shuttle Radar Topography Mission (SRTM) Global with a 30 m raster resolution and downloaded from OpenTopography (https://www.opentopography.org/). The land-use and soil maps from 2003 were obtained from Cambodia National Mekong Committee (CNMC) with a 250-m raster resolution. Weather and hydrological data were collected from the Pursat Provincial Department of Water Resources and Meteorology. There are several rainfall stations nearby and within the river basin, but one station, Kravanh, was selected because it is the only station that has lengthy data recorded from the 1990s and is located in the drained area of the Bak Trakuon outlet (discharge monitoring station used for model evaluation) (Figure 1b). Because of this limitation, a uniform distribution of rainfall in the Kravanh station was applied to the model simulation. For minimum and maximum air temperatures, the dataset in Pursat station located downstream of the river basin was used because it was the only available temperature station in the study area. Other climate data, such as relative humidity, wind speed, and solar radiation, were not available and were simulated by a weather generator. In this study area, the daily river discharge data were the only available hydrological data and were used for model calibration and evaluation.

3.2. SWAT Model Setup

In this study, the ArcSWAT 2012 interface was used to setup the model. With a threshold drainage area of 5000 ha, the basin was discretized into 60 sub-basins, and by classing the slope below and above 5%, the sub-basins were further subdivided into 667 HRUs (Hydrologic Response Units). The model then required the inputs of weather data before it could be run. For this study, the model was run in a monthly time step for a total simulation period of 26 years from 1990 to 2015 in which the period from 1990 to 1994 was treated as the warm-up period to mitigate the initial conditions and was not included in the analysis.

3.3. SWAT-CUP with the SUFI-2 Algorithm

Calibration can be difficult and is almost infeasible for many large-scale applications [7]. A number of auto-calibration and uncertainty analysis tools for SWAT have been developed to solve this problem and are currently available to assist the optimization process. In this study, SWAT-CUP with the SUFI-2 algorithm was used to calibrate the model. SWAT-CUP is an auto-calibration and uncertainty analysis module program developed for SWAT [12,25]. SWAT-CUP allows the use of several algorithms for the calibration or validation procedure. Among them, the SUFI-2 algorithm [16,26], which is a semi-automated approach used for model calibration, validation, sensitivity, and uncertainty analysis, was used in this study. In SUFI-2, all sources of parameter uncertainties are assigned to parameters. The range of input parameters is assumed to be a uniform distribution, while the model output uncertainty is quantified by the 95% prediction uncertainty (95PPU) determined at the 2.5% and 97.5% levels of the cumulative distribution of output variables obtained via Latin hypercube sampling. The details of the SWAT-CUP program and SUFI-2 algorithm are in the SWAT-CUP user manual [12].

3.3.1. Parameterization

Parameterization is a process of selecting parameters and setting their initial value ranges for calibration. Based on previous studies on the region [27,28,29] and literature on SWAT calibration/validation [7,11], nine parameters were selected, and then, their initial value ranges were defined. Two methods, Relative and Replace, were used for different parameters. When Relative is reported this means that a percentage relative change is applied to original parameters, and only when Replace the min and max range values allowed for calibration are reported. The descriptions of each parameter along with their minimum and maximum values are presented in Table 2.

3.3.2. Objective Functions

Eleven objective functions are currently allowed in the SUFI-2 algorithm. Eight of the eleven were selected to evaluate their impacts on the calibration results, calibrated parameters, and the overall discharge processes simulation owing to their popularity and available details on accessibility. For the remaining three objective functions such as the multiplicative form of the square error, the summation form of the square error, and the Chi square, we have no or very little knowledge about them and will not be able to make any discussion and evaluation for their performances. The selected eight objective functions include the coefficient of determination (R²), modified coefficient of determination (bR²), Nash–Sutcliffe efficiency (NSE), modified Nash–Sutcliffe efficiency (MNS), ratio of the standard deviation of observations to the root mean square error (RSR), ranked sum of squares (SSQR), Kling–Gupta efficiency (KGE), and percent bias (PBIAS). The equations and references for these objective functions are presented in Table 3. For NSE and RSE, it is expected that their simulation results are similar as they are equivalent objective functions given that NSE = 1 − RSR² according to their equations in Table 3.

Here, $\mathrm{Q}$ is a variable (e.g., discharge), $\mathrm{m}$ stands for the measured value, $\mathrm{s}$ stands for the simulated value, $\mathrm{i}$ represents the ith measured or simulated data, ${\overline{\mathrm{Q}}}_{\mathrm{m}}$ is the mean measured data of variable $\mathrm{Q}$, ${\overline{\mathrm{Q}}}_{\mathrm{s}}$ is the mean simulated data of variable $\mathrm{Q}$, $\mathrm{b}$ is the coefficient of the regression line between the measured and simulated data, $\mathrm{p}$ is the modified Nash–Sutcliffe efficiency factor ($\mathrm{p}$= 1 was used in this study), $\mathrm{n}$ is the total number of measured or simulated data, $\mathrm{j}$ represents the rank, and $\mathrm{r}$ is the linear regression coefficient between the simulated and measured variables, and $\mathsf{\alpha }=\frac{{\mathsf{\sigma }}_{\mathrm{s}}}{{\mathsf{\sigma }}_{\mathrm{m}}}$ and $\mathsf{\beta }=\frac{{\mathsf{\mu }}_{\mathrm{s}}}{{\mathsf{\mu }}_{\mathrm{m}}}$, where ${\mathsf{\sigma }}_{\mathrm{s}}$ and ${\mathsf{\sigma }}_{\mathrm{m}}$ are the standard deviations of the simulated and measured data, respectively. Additionally, ${\mathsf{\mu }}_{\mathrm{s}}$ and ${\mathsf{\mu }}_{\mathrm{m}}$ are the means of the simulated and measured data, respectively.

3.3.3. Model Calibration, Validation, and Evaluation

Model calibration involves fitting parameter values by comparing the predicted output and measured data until a satisfactory standard for an objective function is achieved [36], whereas validation is an evaluation of whether the calibrated model has the ability to predict for later periods. In this study, the model was run in a monthly time step. The measured discharge data at the Bak Trakuon station during the 1995–2008 and 2009–2015 periods were used for model calibration and validation, respectively.

The strengths of the calibration and validation of each objective function were evaluated. The evaluation was based on statistical indices of all objective functions used in this study. The detailed descriptions and equations are presented in Table 3, and their statistical ranges, optimal values, and satisfactory values are shown in Table 4. The satisfactory threshold values for bR² and MNS were considered greater than or equal to 0.4, which was adopted from Kouchi et al. [20], who introduced measures based on the results of the studies of Muleta [37], Mehdi et al. [38], and Akhavan et al. [39]. For the SSQR, a satisfactory value could not be specified because the measured and simulated variables were independently ranked, and the value depends on the magnitude of the variables being investigated. Satisfactory values of the other indices were based on Moriasi et al. [32] and Thiemig et al. [40].

3.4. Evaluation of Each Objective Function

To identify a better objective function when calibrating the discharge in this study area, the simulated result of each objective function was graphically compared with the observed data, and the performance evaluation criteria (as shown in Table 4) was also taken into consideration. Additionally, characteristics of this tropical monsoon river basin were analyzed, and based on these characteristics, estimated results of water balance components, such as evapotranspiration, surface runoff, lateral flow, groundwater flow, estimated results of the water yield, and calibrated values of the best parameter sets obtained from all objective functions during the calibration process, were also evaluated.

Additionally, the MOD16A2 Version 6 Evapotranspiration/Latent Heat Flux product with an 8-day composite dataset produced at 500 m pixel resolution from MODIS (Moderate Resolution Imaging Spectroradiometer) was used to confirm the actual evapotranspiration amount of the Pursat River Basin [41]. MODIS ET data was obtained for 2001 to 2015 period through R package “MODIStsp” and analyzed using R 4.0.2 version software (The R Foundation for Statistical Computing, Vienna, Austria) [42].

For further evaluation of the hydrological characteristics of this watershed, simulated discharges of the best calibrated parameter sets of the different objective functions were categorized into the base flow, rising limb, peak flow, and falling limb phases and then compared with the observed data at the respective flow phases from 1995 to 2008 (calibration period) using scatter plots. Then, the coefficient of determination (R²) and root mean squared deviation (RMSD) were used for the evaluation. The R², which was calculated as shown in Table 3, determines how much variance the two variables share, and its value varies between 0 (no correlation) and 1 (perfect correlation). However, when the model systematically over- or under-predicts all the time, the R² value is still close to 1. To cope with this problem, the slope and intercept of the regression on which R² is based were taken into account. For a good agreement, the slope and intercept should be close to 1 and 0, respectively. The RMSD represents the mean deviation of the predicted values with respect to the observed values [43,44] and can be calculated as

(1)$\mathrm{RMSD}=\sqrt{\frac{1}{\mathrm{n}-1}{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{Q}}_{\mathrm{m}}-{\mathrm{Q}}_{\mathrm{s}}\right)}_{\mathrm{i}}^2},$

where $\mathrm{n}$ is the total number of observations, ${\mathrm{Q}}_{\mathrm{m}}$ is the measured discharge, ${\mathrm{Q}}_{\mathrm{s}}$ is the simulated discharge, and $\mathrm{i}$ is the ith measured or simulated data. The unit of the RMSD is the same as the unit of ${\mathrm{Q}}_{\mathrm{m}}$ and ${\mathrm{Q}}_{\mathrm{s}}$. The RMSD is always positive, and a value of 0 (almost never achieved in practice) indicates a perfect fit of the data. Generally, a lower RMSD is better than a higher one.

4. Results

4.1. Simulation Results

Figure 3 graphically represents the best simulated discharges on the monthly time step (calibrated from 1995 to 2008 and validated from 2009 to 2015) of different objective functions compared with the observed data. The simulated discharges corresponded well with the rainfall data. With the R² objective function, the calibrated model usually overestimated peak flows. When bR² was used, the peak flow estimation became even more overestimated compared to the observed data. This originated from the equations of these objective functions. Both are defined by a minimization of total errors from a linear regression model and the length from the average value, which are not direct errors of measured and estimated data and could have the effect on the calibrated parameters. According to Legates and McCabe [45], owing to the squared differences in their equations, they are oversensitive to high extreme values and insensitive to additive and proportional differences between model predictions and measured data.

For the choice of NSE, MNS, and RSR objective functions, their calibrated results were similar, especially the NSE and RSR objective functions. In addition to increasing the fitness between the simulated and observed base flows and recession curves, the simulated peak flows were also significantly reduced with considerable correspondence to the observed data. Similarly, when we used the SSQR and KGE objective functions, the simulation results could also capture well the distribution of the observed flow, although the simulated peak flows were slightly high. Those objective functions focus on errors between the measured and estimated data. Then, optimization of these objective functions reduces the errors.

The PBIAS objective function, however, gave a different calibration result. Compared to the calibration results of the other objective functions, the simulated peak flows were occasionally overestimated when we used the PBIAS objective function, whereas the simulated base flows were occasionally underestimated. This is the result of the characteristics of this objective function. When the model overpredicts as much as it underpredicts, PBIAS can still provide a deceptive rating of model performance [35,46].

4.2. Model Performance

The model performances of the calibrated parameter sets obtained from different objective functions were evaluated and compared by using the statistical indices of all objective functions (Figure 4). Many objective functions performed better than satisfactory in terms of calibration and validation for most of the criteria described in Table 4. The simulated results of the fitted parameter sets from the NSE, MNS, RSR, and KGE objective functions satisfied all statistical indices, whereas the results of the parameter set by the R² objective function failed to satisfy the PBIAS statistical index during the validation period. For the result of the parameter set by SSQR, it failed to satisfy the MNS statistical index during both calibration and validation periods. Of the worst, the result of the parameter set by the bR² objective function failed to satisfy NSE, MNS, and RSR statistical indices during the calibration period, whereas the SSQR statistical value of this objective function during the calibration period was quite high compared to that of other objective functions. For the PBIAS objective function, the result of the parameter set during the calibration period failed to satisfy the NSE, MNS, and RSR statistical indices.

In this study, as mentioned earlier, the calibration period was defined longer (1995–2008). The purpose of this was to cover all the magnitudes of discharge pattern which included not only the normal years but also the wettest year in 1996 and the driest year in 2002 (Figure 3) in the calibration process. On the contrary, the validation period, which was shorter from 2009 to 2015, included mostly normal years only. As the result, the performance of the validation as measured by the statistical indices was better than that of the calibration in many cases (Figure 4).

5. Discussion

5.1. Discharge Process Estimations

The results of annual average water balance components and water yields estimated by different objective functions using the best parameter sets in the calibration process are shown in Figure 5. The results showed that the estimated water balance components and water yields differed by the types of objective functions used. The estimated annual average evapotranspiration (ET) by different objective functions ranged from 656.60 to 756.30 mm (48.99% to 56.42% of the total rainfall), whereas the estimated annual average surface runoff (SurQ) ranged from 116.27 to 560.05 mm (8.67% to 41.78% of the total rainfall) (Figure 5a). The ratios of estimated annual average lateral flow (LatQ) and groundwater flow (GWQ) to the total rainfall generated by those objective functions ranged from 0.51% (6.83 mm) to 26.02% (348.76 mm) and 0.31% (4.17 mm) to 13.18% (176.7 mm), respectively. Additionally, the ratios of estimated annual average deep aquifer recharge (DAR) and amount of water moving from the shallow aquifer to plant/soil profile (Revap) to the total rainfall generated by those objective functions ranged from 0.26% (3.51 mm) to 0.89% (11.92 mm) and 0.86% (11.54 mm) to 11.84% (158.66 mm), respectively. For the estimated annual average water yield (Figure 5b), with the MNS objective function, the model generated the lowest annual average water yield of 456.49 mm, whereas the highest annual average water yield of 679.79 mm was given by the model when we used the bR² objective function. The result showed that the estimated water balance components and water yields of the NSE and RSR objective functions were the same.

In a study conducted by Shimizu et al. [47], almost all mountainous regions of west Cambodia (where this study area was also located) had annual renewable freshwater resources (water yield) of 500 mm. Thus, NSE and RSR objective functions provided closer annual average water yield estimations than other objective functions at the same value of 492 mm, followed by the MNS, KGE, and PBIAS objective functions, which generated water yields of 456.49, 544.45, and 548.98 mm, respectively. For the results of the water balance components estimation, of the objective functions used in this study, the MNS, NSE, RSR, PBIAS, and KGE objective functions provided the most reasonable estimation results of annual average evapotranspiration (ET) of 756.3, 726.2, 726.2, 723.9, and 721.7 mm, respectively. These results closely corresponded to the ET value obtained from MODIS ET of 745 mm on average annually (Figure 6) and the findings in the study conducted by JICA [48] in southern Cambodia with an estimated annual average ET of 743.9 mm. The goodness of the ET simulated by these objective functions were also confirmed by the linear regression with the ET obtained from MODIS ET as shown in Figure 7. Again, the PBIAS, NSE, RSR, MNS, and KGE objective functions provided relatively small estimated values of groundwater flow (GWQ), and these objective functions, except for the PBIAS objective function, generated a small proportion of surface runoff (SurQ) compared to a proportion of lateral flow (LatQ), which reflected the flow characteristics of this study area. A detailed description of the discharge characteristics and physiographic condition of the river basin is presented in Section 5.4.

5.2. Best Parameter Sets and Sensitivity Rank

Table 5 shows the fitted parameter values and sensitivity rank of the parameters used in this study obtained from different objective functions during the calibration process. Besides the NSE and RSR objective functions, different objective functions generated different values of best parameter sets and sensitivity rank. The fitted parameter sets of the NSE, RSR, and MNS objective functions were in a reasonable range in this study. When we used the R² objective function, the relative value of the parameter change for CN2 was high (+13%), whereas the existing parameter value of CN2 of the initial model was (between 55 and 92) approximately 78 on average, which was already high for the land-use and soil types of this study area. This led to a big runoff (surface and lateral flow) amount, as shown in Figure 5a, and high simulated peak flows, as shown in Figure 3. This objective function produced a small threshold value of GWQMN of 1680.42 mm, which led to more groundwater flow, and a small coefficient of GW_REVAP of 0.03, which generated a lower evapotranspiration rate and revap because of the limitation of movement of water from the shallow aquifer to the root zone. Moreover, these three parameters were among the five most sensitive parameters, which highly controlled the simulation results of the R² objective function. Even worst, the bR² objective function produced a higher relative value of the parameter change for CN2 of +15%, a smaller threshold value of GWQMN of 291.86, a low coefficient of GW_REVAP of 0.07, and a large value of ESCO, which then led to a large runoff, greater groundwater flow, smaller evapotranspiration, and lower revap (Figure 5a). Consequently, this objective function highly overestimated the peak flows as shown in Figure 3. Furthermore, while the value of the initial model of SOL_AWC was small between 84 and 371 mm (approximately 186 mm on average), the fitted value of the relative change was still underestimated at +16%, and also, the fitted value of SLSUBBSN was overestimated at 68.03 m. As a result, a huge surface runoff and a neglected lateral flow occurred for this objective function. Additionally, the calibrated value of the average slope steepness (HRU_SLP) of bR² was small. However, this parameter was the least sensitive, which may not have as a considerable effect as the earlier mentioned six parameters (they were the top six sensitive parameters).

Regarding the SSQR objective function, the problems were that the fitted value of the relative change of SOL_AWC was small (18%), whereas the fitted value of SLSUBBSN was large (37.25 m), which were the fifth and seventh most sensitive parameters, respectively. This resulted in a large surface runoff and small lateral flow (Figure 5a). Moreover, the large value of ESCO (0.95), which was the third most sensitive parameter, led to a low evapotranspiration; the small value of GW_REVAP (0.04), which was the second most sensitive parameter, led to big groundwater flow and small revap; and the small value of REVAPMN (5.75), which was the sixth most sensitive parameter, contributed to a slightly high deep aquifer recharge. Regarding the KGE objective function, the problems included the calibrated parameter set being slightly large for the relative change of CN2 (8%), which was the fifth most sensitive parameter, the slightly small value of GWQMN (1731.04), which was the third most sensitive parameter, and the slightly small value of GW_REVAP (0.07), which was the most sensitive parameter. As the result, runoff (surface and lateral flow) and groundwater flow generated by this objective function were slightly large, and the revap was relatively small (Figure 5a). For the PBIAS objective function, the calibrated parameter value of the relative change of CN2, which was the most sensitive parameter, was huge (+14%), whereas the value of the relative change of SOL_AWC, which was the second most sensitive parameter, was negative (−10%) and the value of SLSUBBSN, which was the fourth most sensitive parameter, was too long (118.23 m). Consequently, the runoff obtained from this objective function was large with a huge surface runoff and a very small lateral flow (Figure 5a). This is shown in Figure 3 as an overestimation of the peak flows. This objective function also produced a very small value for ESCO, which was the third most sensitive parameter, which should produce a large evaporative demand from the soil. However, owing to the limited available plant water (small SOL_AWC), evaporative demand from the soil was also restricted, leading to an evapotranspiration restriction.

5.3. Hydrograph Components Estimation

To evaluate how different objective functions performed the simulation of each hydrological process, the monthly calibrated results of discharge and observed data were classified into base flow, rising limb, peak flow, and falling limb phases. Based on the monthly average hydrograph in Section 5.4, base flow, which is the period of low flow, was considered from January to March. Then, the rising limb or concentration curve, the ascending portion of the hydrograph, was considered from April to August. Peak flow or crest segment, the inflection point on the rising limb to the falling limb, was considered from September to October; and the falling limb or recession curve, the descending portion from the point of inflection at the end of the crest segment to the base flow, was considered from November to December for each year of the calibration period from 1995 to 2008. In this section, the results of RSR objective function was not further presented and discussed as it is equivalent to and produced the same results with NSE objective function.

Figure 8 presents the scatter plots of the simulated versus observed discharge of base flow period for different objective functions. All of the objective functions generally overestimated the base flows. This behavior can be explained by the small fitted parameter values of GWQMN below 3000 mm as generated by many of the objective functions and the small fitted parameter values of GW_REVAP below 0.05 as defined by some of the objective functions (Table 5). A study conducted by Rafiei Emam et al. [49] in central Vietnam, which is also predominant by forest, defined the final range for GWQMN between 3133 and 3756 mm. However, among them, NSE and MNS objective functions produced better simulation results for the base flows with higher R² values, better slope (optimum value of 1) and intercept (optimum value of 0), and smaller RMSD value. The simulated base flows of the MNS objective function achieved an R² of 0.45, slope of 0.48, intercept of 3.59, and RMSD of 7.78, whereas the simulation results of the NSE objective function provided a value of R² of 0.48, slope of 0.41, intercept of 2.84, and RMSD of 10.76.

For the rising limb estimation, none of the simulation results of these objective functions showed a reasonable correlation with the observed data with an R² of 0.21 at most and large intercepts; however, the correlation slopes of some of them reached a value that was greater than 0.6. Again, NSE objective function provided the closest simulation results (Figure 9) with R², slope, intercept, and RMSD values of 0.21, 0.67, 35.16, and 63.97, respectively.

Conversely, the simulated peak flows of these objective functions gave a reverse performance result (Figure 10). The simulated peak flows of the PBIAS, bR², R², and KGE objective functions attained slightly higher R² values compared to the results of the other objective functions at slightly larger than 0.50, whereas those of the remaining objective functions were a little less than 0.50. However, the regression slopes of the results obtained from the NSE and MNS objective functions reached a satisfactory value of above 0.70, and those obtained from the other objective functions were between 0.54 and 0.64. The intercepts obtained from all objective functions were between 76.82 and 102.54, and the RMSD values were between 106.67 and 127.61.

For the falling limb estimation, the performances of most objective functions were good (Figure 11). The MNS objective function performed well in simulating the falling limbs with an R² of 0.79 with a good slope of 1.01, a small intercept of −9.62, and a RMSD of 37.38. The NSE and KGE objective functions performed similarly but were lower with an R² of approximately 0.78 and RMSD of approximately 41. However, the regression slope and intercept obtained from the NSE objective function were at 0.97 and −14.46, respectively, whereas those obtained from KGE were 0.90 and −5.52, respectively.

The results showed that the NSE and MNS objective functions provided overall better estimation results for all the components of the hydrograph for this river basin. However, KGE, R², bR², SSQR, and PBIAS were among the more poor objective functions, especially for the simulation during low flow periods. For NSE, the differences between the observed and predicted values were calculated as squared values. As a result, larger values in a time series are strongly overestimated, whereas lower values are neglected [45]. Additionally, runoff peaks will tend to be underestimated when NSE is used in the optimization [34]. However, NSE is good for use with continuous long-term simulations and can be used to determine how well a model simulates trends for the output response of concern [46,50]. Because the calibration duration of this study was 14 years, it is likely that the NSE objective function could capture this long-term trend of the discharge. For the MNS objective function, which was the modified form of NSE with the modified factor of p = 1 used in this study, it can be expected that the modified forms are more sensitive to significant over- or under-prediction than the squared forms [30]. However, in this study, the performance of the MNS objective function was only slightly better than the performance of the NSE objective function when we simulated the falling limbs, but it always slightly performed worse than the NSE objective function when we simulated the other components of the hydrograph.

Another objective function used in this study, KGE, is a decomposition of NSE. Similar to NSE, the runoff peaks will tend to be underestimated, but when the KGE optimization is used, the underestimation will not be as severe [34]. As a result, the simulated peak flows of KGE exhibited a slightly better correlation than that of NSE. However, the slope of the regression line for NSE was better (Figure 10), which was because of the suitability of NSE in regressing the observed against the simulated values [34]. The R² objective function is widely used in hydrological modeling studies, but it is oversensitive to high extreme values and insensitive to additive and proportional differences between model predictions and measured data [45]. For the bR² objective function, the under- or over-predictions are quantified together with the dynamics, which results in a more comprehensive reflection of model results [30]. The SSQR objective function aims at fitting the distribution of the flows, ensuring that the full range of the flows is represented but without considering the time of occurrence of a given value of the flows [33]. Perhaps due this characteristic, the estimated recession curves (falling limbs) were typically flatter than the observed data, and the estimated base flows were shortened (Figure 3). As a result, the simulation performances of the falling limbs (Figure 11) and base flows (Figure 8) of this objective function were relatively low. For the PBIAS objective function, it is useful for continuous long-term simulations and can be used to determine how well the model simulates the average magnitudes for the output response of interest [46]. PBIAS can provide a deceptive rating of model performance when the model over-predicts as much as it under-predicts, in which case PBIAS will be close to zero even though the model simulation is poor [46]. This may be why there were several sudden peak and drop points of the simulated result of this objective function during the base flow and rising limb periods (Figure 3), leading to poor model performances in simulating the base flows and rising limbs, as shown in Figure 6 and Figure 7, respectively.

5.4. Objective Functions Corresponding to the Characteristics of the River Basin

Figure 12 shows the monthly average hydrograph of the observed flow at Bak Trakuon Station and rainfall at Kravanh Station during the calibration period from 1995 to 2008 (data between 1997 and 1998 were excluded due to missing data). The validation period was not included because this study focused on the effect of different objective functions on model calibration. Because the dry season, which extends from December to April, is influenced by the northeast monsoon system [22], the river discharge is relatively low, approaching zero between January and March. This indicated that groundwater from the upstream area did not contribute considerably to the river discharge. However the river discharge begins to rise from April when the wet season starts, and the first peak occurs in May as the monsoon rain travels north [22]. This peak is followed by a period of lower rainfall between June and August. The greatest peak occurs between the months of September and October and is caused by a southerly shift in the monsoon circulation pattern, which is characterized by heavy rainfall. This indicated that the flow characteristic in this river basin is highly controlled by the monsoon rainfall pattern, especially during the wettest period from August of November when the correlation between the monthly discharge and the monthly rainfall was so high.

Additionally, the physiographic condition likely influences the hydrological process as well, particularly at beginning of the rainy season from April to July when the correlation between the monthly discharge and the monthly rainfall was low as the result of a lag time in the hydrograph. With the size of the drainage area at Bak Trakuon Station, the elongated shape, the forested land cover with varying densities, and the textures of Dystric Leptosol and Cambisol may also contribute to this broad rising limb in the hydrograph (concentration curve of the hydrograph from April to September) owing to retardation of overland flow and the increase of infiltration and storage capacities of the soils [51]. However, the sharp slope of the falling limb of the hydrograph (recession curve of the hydrograph from late October to December) is likely due to the hilly terrain topography of the drainage area that form a large stream and valley slopes, which result in quick depletion of storage [51].

Compared to other land covers, a forested watershed has some unique features. Mature forests have relatively large aboveground (i.e., over-story and under-story layers) and belowground (i.e., roots) biomass [52,53]. They generally have a higher canopy surface roughness, higher leaf area index, and deeper roots compared to crops and/or grass [54], which result in a relatively high ET [55] and soil infiltration capacity. Forest soil permeability would be maintained by defoliation and organic matter supply from the biomass and considerably reduce the potential for surface flow, lower the total stream flow, and lower the peak flow [56]. Generally, forest stream flow originates from subsurface flow (lateral flow) or groundwater discharge at headwater streams.

Based on the hydrograph and physiographic characteristics of this study area, the NSE, RSR, and MNS objective functions in addition to their goodness of fits for simulating the discharge (Figure 3 and Figure 4) were found to be among the better objective functions for estimating hydrological components, such as a higher ET, lower surface runoff, larger lateral flow, smaller groundwater flow, greater revap, and lower water yield (Figure 5). Moreover, their performances in simulating hydrograph components were overall better than other objective functions, especially the NSE and RSR objective functions, when simulating the base flow, rising limb, and falling limb (Figure 8, Figure 9 and Figure 11).

6. Conclusions

The SWAT model was employed to simulate the stream flow of Pursat River Basin. Eight different objective functions were used in the calibration process of SWAT-CUP with the SUFI-2 algorithm to examine their influences on the calibration results, parameter optimizations, and water resources estimations. Many objective functions performed better than satisfactorily for calibrating the SWAT model. However, different objective functions defined different fitted values and sensitivity ranks of the calibrated parameters, resulting in different estimations of water balance components and water yield. Among the objective functions used in this study, the Nash–Sutcliffe efficiency (NSE) and ratio of standard deviation of the observations to root mean square error (RSR) are equivalent and produced identical simulation results, including the parameter sensitivity and fitted values of the calibrated parameters, leading to the same water balance components and water yields estimation results. By taking results and data of previous studies and according to the characteristics of the river basin, either NSE or RSR objective functions gave the best estimation results of annual average water yield and other water balance components, such as annual average evapotranspiration (ET), groundwater flow (GWQ), surface runoff (SurQ), and lateral flow (LatQ), because these equivalent objective functions generated reasonable fitted values of the calibrated parameters. Moreover, either of them was also better at calibrating the base flow, falling limb, and overall the entire flow phases of the hydrograph for this area.

Author Contributions

Conceptualization and methodological framework, D.S. and T.K.; investigation and data collection, D.S. and P.T.; SWAT model setup, calibration and validation, D.S. and L.H.T.; MODIS ET analysis, A.F. and T.K.; result analysis and visualization, D.S.; writing—original draft preparation, D.S.; writing—review and editing, T.K., C.O., L.H.T., P.T., A.F. and D.S.; supervision, T.K. and C.O. All authors have read and agreed to the published version of the manuscript.

Funding

This work was partially supported by JSPS KAKENHI Grant Number JP18H02295.

Acknowledgments

The authors also would like to thank Enago (www.enago.jp) for the English language review.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. (a) Location of the Pursat River Basin and (b) land-use map in 2003 with monitoring stations.

Figure 2. Workflow diagram.

Figure 3. Monthly simulated and validated discharges of different objective functions compared with observed data.

Figure 3. Monthly simulated and validated discharges of different objective functions compared with observed data.

Figure 3. Monthly simulated and validated discharges of different objective functions compared with observed data.

Figure 4. Results of model performance indices of different objective functions.

Figure 5. (a) Results of the water balance components and (b) water yields of different objective functions using the best parameter sets from the calibration process.

Figure 6. Average monthly (left) and annual evapotranspiration (right) in Pursat River basin (drainage area at Bak Trakuon outlet) from 2001 to 2015 obtained from MODIS ET.

Figure 7. Linear regression between monthly average evapotranspiration obtained from the simulation results (ET Sim.) of different objective functions and monthly average evapotranspiration obtained from MODIS ET.

Figure 8. Scatter plots of simulated versus observed base flows for different objective functions.

Figure 9. Scatter plots of simulated versus observed rising limb for different objective functions.

Figure 10. Scatter plots of simulated versus observed peak flows for different objective functions.

Figure 11. Scatter plots of the simulated versus observed falling limbs for different objective functions.

Figure 12. Monthly average hydrograph at Bak Trakuon Station during the calibration period from 1995 to 2008 (excluding 1997 and 1998 when rainfall data were missing).

Table 1

Description of data used in this study.

Data	Description	Year/Period	Source
Digital Elevation Model (DEM)	Shuttle Radar Topography Mission (SRTM) Global: Raster resolution of 30 m	-	OpenTopography (https://www.opentopography.org/)
Land-use map	Raster resolution of 250 m	2003	Cambodia National Mekong Committee
Soil data	Raster resolution of 250 m	-	Cambodia National Mekong Committee
Weather data	Daily rainfall at Kravanh station	1994–2015	Pursat Provincial Department of Water Resources and Meteorology
	Daily maximum and minimum temperature at Pursat station	2001–2015	Pursat Provincial Department of Water Resources and Meteorology
Hydrological data	Daily discharge at Bak Trakuon station	1995–2015	Pursat Provincial Department of Water Resources and Meteorology

Table 2

Parameters used for the calibration and their initial ranges.

Parameter	Extension	Method	Description	Initial Range
				Min	Max
CN2	.mgt	Relative	SCS runoff curve number	−0.25	0.25
SOL_AWC ()	.sol	Relative	Available water capacity	−0.25	0.25
ESCO	.hru	Replace	Soil evaporation compensation factor	0.01	1
OV_N	.hru	Replace	Manning’s “n” value for overland flow	0.01	30
HRU_SLP	.hru	Replace	Average slope steepness	0	1
SLSUBBSN	.hru	Replace	Average slope length	10	150
GWQMN	.gw	Replace	Threshold depth of water in the shallow aquifer	0	5000
GW_REVAP	.gw	Replace	Groundwater “revap *” coefficient	0.02	0.2
REVAPMN	.gw	Replace	Threshold depth of water in the shallow aquifer for “revap *” to occur	0	500

* Revap is an amount of water moving from shallow aquifer to plants/soil profile in watershed during simulation in mm.

Table 3

List of objective functions used in this study.

Objective Functions	Equation	Reference
Coefficient of determination	${\mathrm{R}}^2=\frac{{\left[{{\displaystyle \sum }}_{\mathrm{i}}\left({\mathrm{Q}}_{\mathrm{m},\mathrm{i}}-{\overline{\mathrm{Q}}}_{\mathrm{m}}\right)\left({\mathrm{Q}}_{\mathrm{s},\mathrm{i}}-{\overline{\mathrm{Q}}}_{\mathrm{s}}\right)\right]}^2}{{{\displaystyle \sum }}_{\mathrm{i}}{\left({\mathrm{Q}}_{\mathrm{m},\mathrm{i}}-{\overline{\mathrm{Q}}}_{\mathrm{m}}\right)}^2{{\displaystyle \sum }}_{\mathrm{i}}{\left({\mathrm{Q}}_{\mathrm{s},\mathrm{i}}-{\overline{\mathrm{Q}}}_{\mathrm{s}}\right)}^2}$	[30]
Modified coefficient of determination	${\mathrm{bR}}^2=\{\begin{array}{c}\left|\mathrm{b}\right|{\mathrm{R}}^2 \mathrm{if} \left|\mathrm{b}\right|\le 1\\ {\left|\mathrm{b}\right|}^{-1}{\mathrm{R}}^2 \mathrm{if} \left|\mathrm{b}\right|>1\end{array}$	[30]
Nash–Sutcliffe efficiency	$\mathrm{NSE}=1-\frac{{{\displaystyle \sum }}_{\mathrm{i}}{\left({\mathrm{Q}}_{\mathrm{m}}-{\mathrm{Q}}_{\mathrm{s}}\right)}_{\mathrm{i}}^2}{{{\displaystyle \sum }}_{\mathrm{i}}{\left({\mathrm{Q}}_{\mathrm{m},\mathrm{i}}-{\overline{\mathrm{Q}}}_{\mathrm{m}}\right)}^2}$	[31]
Modified Nash–Sutcliffe efficiency	$\mathrm{MNS}=1-\frac{{{\displaystyle \sum }}_{\mathrm{i}}{\left|{\mathrm{Q}}_{\mathrm{m}}-{\mathrm{Q}}_{\mathrm{s}}\right|}_{\mathrm{i}}^{\mathrm{p}}}{{{\displaystyle \sum }}_{\mathrm{i}}{\left|{\mathrm{Q}}_{\mathrm{m},\mathrm{i}}-{\overline{\mathrm{Q}}}_{\mathrm{m}}\right|}_{\mathrm{i}}^{\mathrm{p}}}$	[30]
Ratio of the standard deviation of observations to root mean square error	$\mathrm{RSR}=\frac{\sqrt{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{Q}}_{\mathrm{m}}-{\mathrm{Q}}_{\mathrm{s}}\right)}_{\mathrm{i}}^2}}{\sqrt{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{Q}}_{\mathrm{m},\mathrm{i}}-{\overline{\mathrm{Q}}}_{\mathrm{m}}\right)}^2}}$	[32]
Ranked sum of squares	$\mathrm{SSQR}=\frac{1}{\mathrm{n}}{\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}}{\left({\mathrm{Q}}_{\mathrm{j},\mathrm{m}}-{\mathrm{Q}}_{\mathrm{j},\mathrm{s}}\right)}^2$	[33]
Kling–Gupta efficiency	$\mathrm{KGE}=1-\sqrt{{\left(\mathrm{r}-1\right)}^2+{\left(\mathsf{\alpha }-1\right)}^2+{\left(\mathsf{\beta }-1\right)}^2}$	[34]
Percent bias	$\mathrm{PBIAS}=100\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{Q}}_{\mathrm{m}}-{\mathrm{Q}}_{\mathrm{s}}\right)}_{\mathrm{i}}}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{Q}}_{\mathrm{m},\mathrm{i}}}$	[35]

Table 4

Performance evaluations for the monthly discharge simulation [20,32,40].

Indices	R2	bR2	NSE	MNS	RSR	SSQR	KGE	PBIAS
Range	0 to 1	0 to 1	$-\infty$ to 1	$-\infty$ to 1	0 to $\infty$	0 to $\infty$.	$-\infty$ to 1	$-\infty$ to $\infty$
Optimal Value	1	1	1	1	0	0	1	0
Satisfactory Value	>0.5	≥0.4	>0.5	≥0.4	≤0.7	-	≥0.5	<±25

Table 5

Fitted values and sensitivity ranks (values in parentheses) of the calibrated parameters obtained from different objective functions during the calibration process.

Parameter	Fitted Parameter Values and Parameter Sensitivity Ranks (Values in Parentheses) by Different Objective Functions
	R2	bR2	NSE	MNS	RSR	SSQR	KGE	PBIAS
r__CN2.mgt *	13%	15%	3%	−1%	3%	6%	8%	14%
	(4)	(5)	(2)	(2)	(2)	(4)	(5)	(1)
r__SOL_AWC().sol *	45%	16%	45%	37%	45%	18%	38%	−10%
	(6)	(4)	(3)	(5)	(3)	(5)	(6)	(2)
v__ESCO.hru	0.76	0.96	0.76	0.57	0.76	0.95	0.72	0.11
	(3)	(1)	(4)	(4)	(4)	(3)	(4)	(3)
v__OV_N.hru	14.84	10.82	22.38	12.37	22.38	29.61	22.70	27.39
	(9)	(7)	(7)	(9)	(7)	(8)	(8)	(6)
v__HRU_SLP.hru	0.83	0.22	0.74	0.91	0.74	0.96	0.89	0.81
	(7)	(9)	(5)	(3)	(5)	(9)	(9)	(5)
v__SLSUBBSN.hru	13.53	68.03	13.21	11.00	13.21	37.25	11.55	118.23
	(1)	(6)	(1)	(1)	(1)	(7)	(2)	(4)
v__GWQMN.gw	1680.42	291.86	2696.49	2381.61	2696.49	3261.89	1731.04	3006.24
	(2)	(3)	(8)	(6)	(8)	(1)	(3)	(8)
v__GW_REVAP.gw	0.03	0.07	0.11	0.14	0.11	0.04	0.07	0.18
	(5)	(2)	(6)	(8)	(6)	(2)	(1)	(9)
v__REVAPMN.gw	380.62	301.52	133.67	416.76	133.67	5.75	349.10	462.03
	(8)	(8)	(9)	(7)	(9)	(6)	(7)	(7)

* The calibrated CN2 and SOL_AWC values were in relative change, while their original values on average from the initial model were 78 and 186 mm, respectively.