Assessing the sensitivity of SWAT physical parameters to potential evapotranspiration estimation methods over a coastal plain watershed in the southeastern United States

Evapotranspiration (ET) constitutes a key element in developing strategies to optimize the use of water in a watershed and to introduce effective water management practices. It can be expressed as potential ET (PET; the amount of ET if ample water is available) or actual ET (AET). PET is a function of available energy, vapor pressure gradient, and vegetation type (Douglas et al. 2009) while AET is subject to spatial and temporal variability in soil type and moisture, water-table depth, and/or canopy characteristics. PET can be inferred through energy budget methods and/or hydrological water balance (WB) estimation (Samadi et al. 2014), with measurements of vegetation status and climate variables, i.e., stomatal conductance, temperature, solar radiation, air pressure, air humidity, and wind (Mu et al. 2011). PET is first simplified as a reference evapotranspiration (ET₀), which is defined as the evaporation occurring from a land surface with ‘reference crop’ (e.g., green grass of 12 cm; Allen et al. 1998) on soil condition with sufficient water available.

Lu et al. (2005) reported that there are approximately 50 methods for PET estimation, including measurement-based, empirical, mass transfer-based, energy-based, temperature-based, and combination methods. Due to their input data requirements and different assumptions involved in the PET estimation methods, they provide different estimates. Direct measurements, such as those provided by lysimeters, pan evaporation method, eddy flux instrumentation or Bowen-ratio instrumentation, give point values, require constant attendance by skilled personnel, and are based on unverified assumptions (Morton 1983). Empirical PET techniques including temperature-based such as Hargreaves (HG) method, (Hargreaves et al. 1985), radiation-based such as Priestley–Taylor (P–T; Priestley & Taylor 1972) and combination of the two (e.g., Penman–Monteith (P–M); Monteith 1965) provide alternative measures of ET rates, and often offer promising results for the PET losses (e.g., Xu et al. 2006). Both radiation- and temperature-based PET models have been widely used for validating traditional hydrologic models that extrapolated findings at individual sites to the regional scale (e.g., Wilson et al. 2002). At present, empirical techniques have become more practicable in hydrological simulation since PET is a common input for watershed models that offers an upper limit to evapotranspirative water losses. Although a clear definition of the ‘best’ method for hydrologic computation is not still evident and the method choice is often subjective (e.g., Xu & Singh 2002; Douglas et al. 2009), primarily depending on the objectives of the study and the type of data available (Verstraeten et al. 2008).

A large body of literature related to the application of PET estimation methods exists in hydrology communities. For instance, Almorox & Grieser (2016) compared the Hargreaves–Samani (HS) method with the P–M equation (FAO-56) and calibrated the coefficients of HS for different climates as represented by the Köppen classification. Their results suggest that the correlation between long-term monthly means of HS and FAO-56 can be improved significantly by using climate-class specific coefficients. In another study, Koedyk & Kingston (2016) investigated the impact of error and uncertainty associated with PET on river flows using climate change scenario in the Waikaia River, New Zealand. They found less sensitivity of runoff projections to PET methods that needs to be explored more in the future. Alemayehu et al. (2015) studied the use of different sources of weather input data and PET estimation methods to derive ET in Kenya. They found that WB components such as ET, deep aquifer loss, and reevaporation from the shallow aquifer vary in magnitude, depending on the data and methods being used. Their study also showed that the HG method is more robust and realistic compare to the P–T and the P–M models.

Douglas et al. (2009) evaluated observed daily ET at 18 sites in Florida (USA) and stressed that the P–T performance is superior to the other methods (P–M and the Turc method (Tc)) for a variety of land covers. Weiß & Menzel (2008) compared the P–T method with the P–M and the HG methods, and advocated that the P–T results were closest to available pan evaporation data. Oudin et al. (2005) studied a total of 23 methods for PET estimation using a variety of micrometeorological input data. They compared the impact of the PET selection methods on four rainfall–runoff models over 308 watersheds located in France, Australia, and the United States. Their results showed that temperature-based PET estimates perform as well as or better than more physically based PET methods. Similar to Oudin et al. (2005), Vörösmarty et al. (1998) compared the performance of 11 different PET methods (simple to more physically based approaches) and achieved similar results. However, Oudin et al.'s (2005) study focused on systematic bias reduction using the P–M method prior to use in the rainfall–runoff models. They used relative performance of the PET methods in rainfall–runoff modeling and compared the computed PET with empirical ET measurements as subjective approaches.

This study focuses on the relative performance of the PET estimation methods in a distributed rainfall–runoff modeling and addresses the impact of the three different PET methods (i.e., P–T, HG, and P–M) on the watershed physical parameters over a forested wetland dominated ecosystem in the southeastern United States (SEUS). A few studies have used hydrology models to compare the relative performance of the PET selection methods in SEUS. For instance, Lu et al. (2005) compared several PET methods across 36 forested watersheds in SEUS with the result that the P–T method performed better than the rest of the approaches. Rao et al. (2011) explored different ways to identify appropriate PET models for two small forested watersheds in SEUS. They stressed that different PET methods provided significantly different estimates, but P–T resulted in the most reasonable estimates of forest ET rates.

While these studies provide insight into the PET estimation method, a challenge to selecting a representative PET method for heterogeneous regions is how the sensitivity of physical parameters may transfer/change depending on the method being used, and how (quantity) the sensitivity of physical parameters may affect the simulated flow regime. Due to scaling error in the PET estimation (e.g., Williams et al. 2004; Ford et al. 2007; Ngongondo et al. 2013) and often too simply derived as the residual of the WB (Domec et al. 2012), PET estimation methods have a large control on the sensitivity of physical parameters in rainfall–runoff modeling. However, although the impact of different PET estimation methods on the watershed physical parameters has been explored before (Andréassian et al. 2004; Zhao et al. 2013; Tamm et al. 2016; Singh et al. 2016; Koedyk & Kingston 2016), few studies have focused on the sensitivity of physical parameters to the PET estimation methods using uncertainty algorithm (e.g., Samadi et al. 2014; Alemayehu et al. 2015).

The objective of this study was to characterize the sensitivity of the Soil and Water Assessment Tool (SWAT) to the PET estimation methods using Bayesian uncertainty algorithm. The approach was to calibrate the three SWAT PET models (i.e., SWAT(P–T), SWAT(P–M), SWAT(HG)) using sequential uncertainty fitting (SUFI2) framework and assess the sensitivity of SWAT physical parameters during a simulation period over a heterogeneous watershed system in the SEUS. The results are compared in terms of their ability to capture the magnitudes of observed streamflow and flow duration curves (FDC). A unique aspect of this research is that this study is the first attempt to assess the sensitivity of physical parameters to the selected PET estimation methods in the SEUS.

The next section presents a description of the study area and introduces the SWAT model, the SUFI2 algorithm, and the PET methods. The methodology is then applied in the following section to identify sensitive parameters and to calibrate streamflow of a distributed rainfall–runoff model. The final section discusses and summarizes the findings of this research.

This watershed typically experiences high climate variability and hydrological changes especially during wet and dry intervals. The climate of the watershed is humid subtropical in which winters are normally mild and summers are long and hot. Precipitation in the summer is dominated by convection storms, and in the winter by frontal boundaries. The long-term annual average temperature (1946–2009) is 16.88 °C, and averaged annual precipitation was a little more than 1,300 mm in the same period. Maximum temperature during the same period was 41.66 °C while minimum temperature was −14 °C.

The major soil type of this landscape is sandy loam with moderate permeability (see soil map in Figure 1). According to the SWAT output, the proportion of low storage shallow soil is around 90% which is restricted in moderately well drained (hydrologic group B soil) to poorly drained soil (hydrologic group D soil) with agriculture, rangeland, forested rangeland, and forested wetland land covers. Table 1 exhibits the percentage for each land use class (land use data of 2006) in the Waccamaw watershed. Forested wetland and evergreen forest are two major land use land cover classes in the study area.

Sandy and sandy loam soils are the dominant types in the watershed. Approximately 67% of the soils are one of the four series, all of which are either hydrologic group B, D, or B/D. Only 7% of the soil is hydrologic group A; there is no group C soil (Table 2). Hydrologic group D soil (poorly drained) is adjacent to the main channel and hydrologic group B and B/D soils (mostly of the Rains series in the downstream reaches) are adjacent to those.

Model calibration, validation, and sensitivity analysis was done using data from US Geological Service (USGS) gauging station at Longs. Research data were retrieved on a daily basis from the National Climatic Data Center and USGS portals on 6 September and 13 September 2012, respectively. Selected climate stations (Loris, Whiteville 7, and Longwood) are distributed evenly throughout the watershed.

According to Monteith (1965), during ET processes water vapor is transported at a rate equal to the product of the water vapor content and the wind speed across a landscape. As shown in Figure 2, when air moves from a dried surface (upland) to a wet area (riparian buffer) in the coastal plain watershed, the concentration of water vapor increases at the transition to a higher value downwind. Dense vegetation and rich canopy as well as saturated/semi-saturated soil help to raise the concentration values throughout alluvial riparian floodplains. At the transition level, the AET rate may immediately maximize to a much higher value compared to that over the dry land, and then it may decline slowly to a value representative of the wetter condition. The lower AET over the dry surface in summer means the overpassing air may be hotter and drier, thus increasing the available heat energy to increase evaporation in the downwind wetter area (e.g., Morton 1983). These processes are a simple illustration of the AET processes in a coastal plain system and the context of this process is much more complicated in reality.

SWAT offers three options for estimating PET using well-known empirical models, i.e., HG, P–T, and P–M. Users also can incorporate measured ET values in the model. A new CN method that includes a one-parameter soil moisture depletion curve to account for antecedent soil moisture conditions was used in this study to compute surface runoff. SWAT computes evaporation as a function of PET by combining ET values from soil and plants separately (e.g., Ritchie 1972). SWAT first evaporates any rainfall intercepted by the plant canopy. Next, it calculates the maximum amount of transpiration and sublimation/soil evaporation. The actual amount of sublimation and evaporation from the soil surface is then calculated. If snow is presented in the HRU, sublimation occurs and when there is no snow (such as this case study), only evaporation from the soil surface is calculated (Neitsch et al. 2005). The same procedure will be demonstrated for canopy and plant transpiration. Therefore, ET is the sum of the evaporation from the canopy and the soil as well as plant transpiration (Alemayehu et al. 2015).

This research examined a Bayesian framework (i.e., SUFI2) for the SWAT parameter sensitivity analysis. SUFI2 uses a multiple regression analysis to compute the statistics of parameter sensitivity. SUFI2 sensitivity analysis defines which parameters have significant influences on the model simulations by presenting the smallest P-value (values close to zero suggest high level of significance) and the largest t-stat (a measure of sensitivity, larger in absolute values are more sensitive). The t-stat represents the precision of the regression coefficient. It is the coefficient of a parameter divided by its standard error. If a coefficient is large compared to its standard error, then it is probably different from 0 and the parameter is sensitive. SUFI2 then compares the t-stat of a parameter with the values in the Student's t-distribution (defines the behavior of the mean of a sample) table to determine the P-value. The P-value for each term tests the null hypothesis that the coefficient is equal to zero (no effect). A low P-value (<0.05) indicates that the null hypothesis should be rejected. Thus, parameters with low P-values are likely to be meaningful additions to the model. Conversely, parameters with large P-values are more weakly associated with changes in the response (the parameter is not sensitive). Table 3 depicts the 15 parameters selected in the model calibration process for the SWAT model.

SUFI2 first defines the goal function with meaningful parameter ranges ⁠. Then, a Latin hypercube sampling is carried out in the initial hypercube to evaluate the corresponding goal function and to perform the sensitivity analysis (based on sensitivity matrix and the parameter covariance matrix). The 95PPU (calculated at the 2.5% and 97.5% levels of the cumulative distribution of output variables) of a sensitive parameter is then computed to find the best solutions (i.e., parameters which produce the optimal goal function). The 95PPU and model evaluation indices are then calculated. If model performances have satisfactory values, the parameter range travels from to as the posterior parameter distribution. Based on Bayes’ theorem, the probability density of the posterior parameter distribution will be driven from the prior density function to calculate the likelihood function of the model (see Yang et al. 2008). Uncertainty is then defined as discrepancy between measured and modeled variables where it is counted by variation between them.

Each SWAT model was calibrated using SUFI2 against five years of daily streamflow data (2003–2007), excluding the three-year warming-up period (i.e., 2000–2002). The calibration of each model was set for four iterations (each iteration representing 500 simulation runs). In each iteration, prior parameter ranges were updated by calculating the sensitivity matrix and the equivalent of a Hessian matrix, followed by the calculation matrix. Parameters were then updated in such a way that new ranges were always smaller than previous ranges and were centered on the best simulation.

In this study, the basic statistics used in the evaluation of the model performance are: Kling–Gupta efficiency (KGE; Gupta et al. 2009), Nash–Sutcliffe efficiency (NSE; Nash & Sutcliffe 1970), coefficient of determination R² multiplied by the coefficient of the regression line (BR²), and sum of the squares of the differences of the measured and simulated values after ranking (SSQR).

The SSQR is highly influenced by the magnitudes of the differences between the observed and simulated values. n is the number of pairs and i represents the rank. Overall, SSQR assesses the relative magnitudes of the bias in the model.

In this study, global sensitivity analysis was implemented for SWAT sensitivity test. This was conducted by allowing all parameters to be changed from their absolute values during the simulation period. This research evaluated and compared parameter sensitivity using a five-year calibration period with different starting points. Each simulation period was shifted by one year, such that subsequent periods have four years of data in common. Overall, ten different calibration periods were considered, and for each dataset, parameters were inferred using the SUFI2 model. This experience suggests that five years of daily streamflow data contain enough information about the parameters estimation of SWAT inputs using different PET estimation methods. Therefore, no significant variation in parameter sensitivity exists between different calibration periods for each PET setup model.

The sensitivity analysis revealed that when the PET estimation method changes, the sensitivity ranking of the parameters varies significantly. P–M method showed more sensitivity to ALPHA_BF, SLSUBBSN, and CH_K2 whereas HG is most sensitive to ESCO, CH_N2, and CN. Sensitivity results of the P–T method further demonstrated that the SWAT model is highly sensitive to SOL_K, CH_N2, and CH_K2 (see Table 4). Differences in sensitive parameters imply that the impact of the PET estimation methods on the ranking of the SWAT sensitive parameters is relatively high for the study watershed. In other words, changing the PET method has a great effect on physical parameters as well as streamflow simulation in this watershed. This could be attributed to the fact that the magnitudes and temporal dynamics of the streamflow showed a high sensitivity to physical parameters controlled by the PET values. Given the differences in the magnitudes of the estimated streamflow, SWAT parameters such as SOL_K, ESCO, and CN show a remarkable variation in rank by changing the PET estimation method. Because these parameters control physical processes and surficial fluxes in the watershed, these results stress that appropriate selection of the PET estimation method is important for appropriate representation of physical processes as well as estimating streamflow and WB components in the watershed under study.

In addition, the error in streamflow prediction (SSQR; not shown in Figure 5) using different PET methods ranged from 61.43 (cubic meter per second; cm/s) to 271.55 cm/s and 294.57 cm/s, respectively, for P–T, HG, and P–M. HG and P–M methods yielded the highest error compared to the P–T method. This may be partly attributed to overpredicted/underpredicted ET during 2007 when the moisture was limited. The large discrepancy in 2007 was related to rainfall deficiency (>− 480 mm below average) and drought condition during the water year (1 October 2006 to 30 September 2007). However, spatial errors in rainfall data during the simulation period may also increase error and bias. For example, the average rainfall measured at the Loris station in 2006 was 1,253.8 mm compared to 1,580.3 mm at the Longwood gauge.

As expected, May to August is the peak time of ET in the Waccamaw watershed. It is interesting to note that AET was highly maximized during May to August 2005 with July having the highest AET rate (>277 mm; Figure 8(a)). Due to wet antecedent soil condition during 2003 (above normal rain of 1,580 mm) and somewhat during 2004 (normal rain of 1,240 mm but with wet antecedent condition), annual AET of 2005 was the lowest value in the study area. The annual AET remained close to the PET rate (around 90% of PET) for the years exceeding the long-term average rainfall and/or the years with just below the average but with the wet antecedent year. It appears in most of the year, the P–T method calculated moderate AET/PET values. In contrast, HG overestimated AET compared to P–T and P–M during summer season while underestimating during winter. The P–M method, on the other hand, overestimated AET during winter and underestimated during summer. From October 2006 to the end of 2007 with consistently below average monthly rainfall, all methods yielded high monthly AET with July 2007 being the highest value (>265 mm). A similar result was obtained by Lu et al. (2005), Douglas et al. (2009), and Rao et al. (2011). Wang et al. (2006) stressed that the sensitivity of the simulated AET to the choice of the PET estimation method were insignificant at a 5% significance level for the southeastern USA. Furthermore, Amatya & Jha (2011) advocated that SWAT is more sensitive to the choice of the PET estimation method; hence, it is important to select the most accurate and spatially stable model for the coastal plain watershed.

As shown in Figure 8(b) and 8(c), the magnitudes of the simulated potential plant transpiration and potential soil evaporation depend on the choice of the PET estimation method. Although potential plant transpiration estimation method is different in the PM, HG, and P–T methods, P–M estimated higher potential plant transpiration compared to the rest of the methods. Nevertheless, relatively lower is obtained from the simplified and less data-dependent method, i.e., the HG model. Overall, the annual average actual soil and plant transpiration reflects the effect of the choice of the PET estimation method on these physical properties, whereas the P–M method results in higher estimates in most cases (Figure 8(b)–8(d)). On average, as shown in Figure 8(d), the P–T method demonstrated a higher annual evaporation from canopy storage, despite the relatively moderate PET estimates, as explained earlier. This is attributed to a higher EPCO value in the P–T method, as shown in Figure 4. The results for the EPCO parameter reveal the sensitivity of the simulated plant uptake compensation factor to the choice of the PET estimation method. Higher EPCO reveals that most of the transpiration demands are met by lower layers in the soil.

The field capacity is the amount of water remaining in the soil a few days after having been wetted and after free drainage has ceased. It appears this value is overestimated by the P–M method while HG underestimated it compared to P–T. The total available water (holding) capacity is the proportion of water available, stored, and/or released between field capacity and the wilting point water contents. The soil types with higher total available water content are generally more conducive to high biomass productivity because they can provide adequate moisture to plants during water shortage. Soils with low available water capacity may need more irrigation and this is the case for shallow rooted crops (i.e., potatoes) located on sandy soil in the study area. Higher available water capacity leads to increased soil moisture in the soil profile, which in turn reduces the initial abstractions, hence, higher surface runoff generation.

In addition, analysis revealed that the highest ET rate is calculated for subwatersheds 25 and 5 with forested wetland land cover due to its potential for ET processes while the least ET is calibrated, respectively, for range and agricultural lands at subwatersheds 26 and 19. Spatial distribution of ET estimation methods also indicates that the HG model overestimated ET for forested wetland, agricultural, and rangeland land covers compared to the rest of the methods (for instance, see HG estimated values for subwatersheds 20 and 25). All three methods also calculated the same values for the forest land cover. Overall, ET distribution showed that agricultural land covers presented moderately higher ET values. Furthermore, all methods presented approximately similar values for groundwater contribution at each subwatershed level.

The aim of this study was to examine the sensitivity of SWAT physical parameters to changes in the PET methods over a coastal plain ecosystem in the southeastern USA. The output of SWAT was compared in terms of parameter sensitivity and streamflow simulations using Bayesian uncertainty algorithm.

The use of ET equations in the SWAT model cause different degrees of sensitivity (and consequently increasing bias) for physical parameters because of either the assumptions made in interrelationship between climatic variables associated with the PET equation or simplified relationship between atmospheric variables, soil and plant/crop computations. Analysis shows that among all three PET estimation models examined here, radiation-based P–T equation which was developed for warm and humid climate condition, performed best for the Waccamaw coastal plain watershed when judged by streamflow simulation, sensitivity analysis, and physical parameter estimation. This reveals that the radiation-based method is efficient (solar radiation, air temperature, and relative humidity control ET process); hence the P-T model may be preferable to the HG and P–M models for the coastal plain with wet and humid surface, as stressed by Lu et al. (2005), Wang et al. (2006), and Rao et al. (2011).

Although there are large differences in PET estimates by various methods during the dormant season (winter) compared to the growing season (spring and summer), modeling results revealed insignificant differences among these three methods in respect to AET estimates. Further, PET estimates remarkably differed during peak and low flow events as well as winter and summer seasons. This analysis reveals that low and high flows are more sensitive to the changes in the PET estimation method. Thus, a careful selection of the PET method is required when using SWAT for extreme simulation (drought and flood) in the study watershed.

The adequacy of the assumptions made in the P–T equation for the Waccamaw watershed has been validated by a number of independent investigations before (Lu et al. 2005; Wang et al. 2006; Samadi et al. 2014) in which it was commonly found that, in vegetated areas with no or very small water deficiency, approximately 95% of the annual evaporative demand is supplied by radiation (Stagnitti et al. 1989). This research finding is also consistent with the study by Rao et al. (2011), who compared three common PET models (FAO-56 grass reference ET, Hamon PET, and P–T) in the southeastern Appalachian mountain region. They found that the radiation-based equation (i.e., P–T method) generally performed better than those that included only temperature-related input variables.

Overall, this results of this study showed that AET represents 35%–53% of precipitation losses in hydrologic modeling of the Waccamaw coastal plain watershed. It is fairly in agreement with Palmroth et al. (2010), who stressed that 30% of AET was explained by the variation in precipitation in the coastal plain watershed. However, the quality of input data and methods for calculating PET (and hence, AET) matters and should receive more attention when applying hydrology models in data-scarce coastal plain watersheds. It is recommended that a comparison of SWAT and remotely sensed ET estimates using climate change scenarios (e.g., Awan et al. 2016; Tamm et al. 2016) be made to provide spatial distributed ET values at subwatershed or HRU level.

The author appreciates those people and agencies that assisted in accessing research data and information. This research was partially supported by South Carolina Sea Grant Consortium (Grant 15520-GA11) and the University of South Carolina (Grant 15520-16-40787). I also wish to thank the anonymous reviewers and the editor for their valuable comments.