One of the key inputs of a hydrologic budget is the potential evapotranspiration (PET), which represents the hypothetical upper limit to evapotranspirative water losses. However, different mathematical formulas proposed for defining PET often produce inconsistent results and challenge hydrological estimation. The objective of this study is to investigate the effects of the Priestley–Taylor (P–T), Hargreaves, and Penman–Monteith methods on daily streamflow simulation using the Soil and Water Assessment Tool (SWAT) for the southeastern United States. PET models are compared in terms of their sensitivity to the SWAT parameters and their ability to simulate daily streamflow over a five-year simulation period. The SWAT model forced by these three PET methods and by gauged climatic dataset showed more deficiency during low and peak flow estimates. Sensitive parameters vary in magnitudes with more skew and bias in saturated soil hydraulic conductivity and shallow aquifer properties. The results indicated that streamflow simulation using the P–T method performed well especially during extreme events’ simulation.

## INTRODUCTION

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 _{0}*), 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.

## METHODOLOGY

### Study area

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.

Area | |||
---|---|---|---|

Land use type | Land use description | sqkm | % |

WATR | Open water | 48.71 | 1.6 |

URLD | Rural residential | 112.81 | 3.7 |

URMD | Medium density residential | 28.68 | 0.9 |

URHD | High density residential | 4.40 | 0.1 |

SWRN | South western range | 1.05 | 0.1 |

FRSD | Forest and woodland | 4.21 | 0.1 |

FRSE | Evergreen forest | 2.48 | 27.0 |

FRST | Mixed forest | 843.97 | 1.4 |

RNGB | Range shrub land | 42.20 | 10.8 |

RNGE | Grasslands/Herbaceous | 334.82 | 5.0 |

HAY | Pasture | 154.90 | 0.1 |

AGRR | Agricultural land-row crops | 2.49 | 20.5 |

WETF | Forested wetland | 641.75 | 27.7 |

WETN | Non-forested wetlands | 864.11 | 1.0 |

UIDU | Industrial | 28.99 | 0.1 |

Area | |||
---|---|---|---|

Land use type | Land use description | sqkm | % |

WATR | Open water | 48.71 | 1.6 |

URLD | Rural residential | 112.81 | 3.7 |

URMD | Medium density residential | 28.68 | 0.9 |

URHD | High density residential | 4.40 | 0.1 |

SWRN | South western range | 1.05 | 0.1 |

FRSD | Forest and woodland | 4.21 | 0.1 |

FRSE | Evergreen forest | 2.48 | 27.0 |

FRST | Mixed forest | 843.97 | 1.4 |

RNGB | Range shrub land | 42.20 | 10.8 |

RNGE | Grasslands/Herbaceous | 334.82 | 5.0 |

HAY | Pasture | 154.90 | 0.1 |

AGRR | Agricultural land-row crops | 2.49 | 20.5 |

WETF | Forested wetland | 641.75 | 27.7 |

WETN | Non-forested wetlands | 864.11 | 1.0 |

UIDU | Industrial | 28.99 | 0.1 |

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.

Hydrologic | Area | |||
---|---|---|---|---|

Series name | soil group | Drainage class | sqkm | % |

Meggett | D | Poorly drained | 559.854 | 18.0 |

Bohicket | D | Very poorly drained | 0.004 | 0.0 |

Kenansville | A | Well drained | 49.424 | 1.6 |

Leaf | D | Poorly drained | 144.951 | 4.7 |

Woodington | B/D | Poorly drained | 339.478 | 11.1 |

Croatan | D | Very poorly drained | 473.114 | 15.2 |

Rains | B/D | Poorly drained | 614.282 | 19.7 |

Norfolk | B | Well drained | 447.622 | 14.4 |

Autryville | A | Well drained | 89.034 | 2.9 |

Leon | B/D | Poorly drained | 42.46 | 1.5 |

Kureb | A | Excessively well drained | 15.614 | 0.5 |

Trebloc | D | Poorly drained | 163.43 | 5.2 |

Nansemond | A | Moderately well drained | 23.525 | 0.7 |

Montevallo | D | Well drained | 35.07 | 1.1 |

Osier | D | Poorly drained | 0.26 | 0.0 |

Chisolm | A | Well drained | 118.71 | 3.8 |

Hydrologic | Area | |||
---|---|---|---|---|

Series name | soil group | Drainage class | sqkm | % |

Meggett | D | Poorly drained | 559.854 | 18.0 |

Bohicket | D | Very poorly drained | 0.004 | 0.0 |

Kenansville | A | Well drained | 49.424 | 1.6 |

Leaf | D | Poorly drained | 144.951 | 4.7 |

Woodington | B/D | Poorly drained | 339.478 | 11.1 |

Croatan | D | Very poorly drained | 473.114 | 15.2 |

Rains | B/D | Poorly drained | 614.282 | 19.7 |

Norfolk | B | Well drained | 447.622 | 14.4 |

Autryville | A | Well drained | 89.034 | 2.9 |

Leon | B/D | Poorly drained | 42.46 | 1.5 |

Kureb | A | Excessively well drained | 15.614 | 0.5 |

Trebloc | D | Poorly drained | 163.43 | 5.2 |

Nansemond | A | Moderately well drained | 23.525 | 0.7 |

Montevallo | D | Well drained | 35.07 | 1.1 |

Osier | D | Poorly drained | 0.26 | 0.0 |

Chisolm | A | Well drained | 118.71 | 3.8 |

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.

### ET processes in the coastal plain 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 model

*et al.*(1993). The SWAT system is embedded within a geographic information system that can integrate various spatial environmental data including soil, land cover, climate, and topographic features. Based on this spatial information, a watershed is discretized into subwatersheds, and these in turn are subdivided into hydrological response units (HRUs). In a distributed hydrology model, a HRU is a unique area that is created by combination of land use, soil, and slope classes. These three classes have large effects on physical properties and runoff simulation in the model. Each process is treated at the HRU level on a

*t*time step to provide aggregated inputs into a river network (Neitsch

*et al.*2011). SWAT then computes five different storages values at each HRU including snow, canopy storage, the soil profile with up to ten layers, a shallow aquifer, and a deep aquifer. SWAT is based on the WB equation (Equation (1)): where Δ

*S*is changes in water storage;

*t*denotes daily time step; and

*P*is the daily amounts of precipitation.

*Q*denotes the total amount of water yield.

_{total}*ET*is evapotranspiration and

*Losses*represents water losses via deep aquifer, percolation, and the evaporation from the shallow aquifer (REVAP).

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).

### PET estimation methods in the SWAT model

*et al.*2015). This method is usually used to estimate potential evaporation from a vegetated surface. PET can be estimated by the P–M method using the following formula (Equation (2)): where

*PET*represents potential evapotranspiration from a reference crop (mm/day); Δ is the slope of saturation vapor pressure–temperature curve (kPa/°C); denotes net radiation (MJ/m

^{2}/day);

*G*is heat flux density to the ground (MJ/m

^{2}/day); is air density (kg/m

^{3}); denotes specific heat at constant pressure (MJ/kg/°C); represents saturation vapor pressure of air at height z (kPa); is water vapor pressure of air at height z (kPa); γ denotes psychrometric constant (kPa/°C); is plant canopy resistance or the (bulk) surface resistance (s/m); and is diffusion resistance of air layer (aerodynamic resistance; s/m).

*PET*). The mass transfer term in the Penman combination equation approaches zero and the radiation terms dominate. Priestley & Taylor (1972) found that the AET from well-watered vegetation was generally higher than the equilibrium potential rate and could be estimated by multiplying the

_{eq}*PET*by a factor equal to 1.26: where λ denotes the latent heat of vaporization (MJ/kg),

_{eq}*E*represents the potential evapotranspiration (mm d

_{o}^{−1}), is a coefficient, Δ indicates the slope of the saturation vapor pressure–temperature curve, de/dT (kPa/°C), γ signifies the psychrometric constant (kPa/°C), is the net radiation (MJ/m

^{2}/d), and

*G*represents the heat flux density to the ground (MJ/m

^{2}/d). Although, the value of may vary throughout the day (Munro 1979), there is general agreement that a daily average value of 1.26 is applicable in a humid climate (Stewart & Rouse 1976; De Bruin & Keijman 1979; Shuttleworth & Calder 1979). Generally, the coefficient for an expansive saturated surface (like the coastal plain watershed) is usually greater than 1.0. This means that true equilibrium PET rarely occurs; there is always some component of advection energy that increases the actual evapotranspiration. Higher values of , ranging up to 1.74, have been recommended for estimating PET in more arid regions (ASCE 1990).

*et al.*1985). This method has a minimum weather data requirement, using only the maximum and minimum air temperature. The equation claims to self-compensate for the lack of solar radiation, humidity, and wind run data and to implicitly account for the effects of cloudiness, relative humidity, and vapor pressure deficit (Hargreaves & Allen 2003). The HG method equation is formulated as follows (Equation (4)): where is incoming extraterrestrial solar radiation (MJ/m

^{2}/d). and are respectively, maximum and minimum air temperature for a given day (°C). When the HG method is used in SWAT, an adjusted evaporative water demand is first calculated by subtracting the free water evaporation from the canopy storage :

*E*) is estimated as a function of the simulated leaf area index (LAI): The potential soil evaporation is a function of and the soil cover index : The soil cover index depends on the above-ground biomass and residue (kg/ha). The actual soil evaporation is limited by the soil water content and is reduced exponentially when the soil water content drops below field capacity.

_{t}### Model calibration and sensitivity analysis

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 signiﬁcant 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.

SWAT parameter^{a} | Name of SWAT parameter | Absolute values | Parameter type | Identifier code |
---|---|---|---|---|

CH_N2 | Manning's ‘n’ value for the main channel | −0.01–0.3 | Channel route (.rte) | Replace^{b} |

SOL_K | Saturated hydraulic conductivity (mm/hr) | 0–2,000 | Soil (.sol) | Relative |

CH_K2 | Effective hydraulic conductivity in tributary channel alluvium (mm/hr) | −0.01–500 | Channel route (.rte) | Replace |

SLSUBBSN | Average slope length (m) | 10–150 | HRU (.hru) | Replace |

ESCO | Soil evaporation compensation factor | 0–1 | HRU (.hru) | Replace |

OV_N | Manning's ‘n’ value for overland flow | 0.01–30 | HRU (.hru) | Replace |

SOL_BD | Moist bulk density (Mg/m^{3} or g/cm^{3}) | 0.9–2.5 | Soil (.sol) | Relative |

GWHT | Initial groundwater height (m) | 0–25 | Groundwater (.gw) | Relative |

EPCO | Plant uptake compensation factor | 0–1 | HRU (.hru) | Replace |

SHALLST | Initial depth of water in the shallow aquifer (mm H_{2}O) | 1,000–50,000 | Groundwater (.gw) | Relative |

ALPHA_BF | Base flow alpha factor (days) | 0–1 | Groundwater (.gw) | Replace |

Curve number (CN) | SCS runoff CN for moisture condition II | 35–98 | Groundwater (.gw) | Relative |

RCHRG_DP | Deep aquifer percolation fraction | 0–1 | Groundwater (.gw) | Replace |

GW_REVAP | Groundwater ‘revap’ coefficient | 0.02–0.2 | Groundwater (.gw) | Replace |

GW_DELAY | Groundwater delay time (days) | 0–500 | Groundwater (.gw) | Replace |

SWAT parameter^{a} | Name of SWAT parameter | Absolute values | Parameter type | Identifier code |
---|---|---|---|---|

CH_N2 | Manning's ‘n’ value for the main channel | −0.01–0.3 | Channel route (.rte) | Replace^{b} |

SOL_K | Saturated hydraulic conductivity (mm/hr) | 0–2,000 | Soil (.sol) | Relative |

CH_K2 | Effective hydraulic conductivity in tributary channel alluvium (mm/hr) | −0.01–500 | Channel route (.rte) | Replace |

SLSUBBSN | Average slope length (m) | 10–150 | HRU (.hru) | Replace |

ESCO | Soil evaporation compensation factor | 0–1 | HRU (.hru) | Replace |

OV_N | Manning's ‘n’ value for overland flow | 0.01–30 | HRU (.hru) | Replace |

SOL_BD | Moist bulk density (Mg/m^{3} or g/cm^{3}) | 0.9–2.5 | Soil (.sol) | Relative |

GWHT | Initial groundwater height (m) | 0–25 | Groundwater (.gw) | Relative |

EPCO | Plant uptake compensation factor | 0–1 | HRU (.hru) | Replace |

SHALLST | Initial depth of water in the shallow aquifer (mm H_{2}O) | 1,000–50,000 | Groundwater (.gw) | Relative |

ALPHA_BF | Base flow alpha factor (days) | 0–1 | Groundwater (.gw) | Replace |

Curve number (CN) | SCS runoff CN for moisture condition II | 35–98 | Groundwater (.gw) | Relative |

RCHRG_DP | Deep aquifer percolation fraction | 0–1 | Groundwater (.gw) | Replace |

GW_REVAP | Groundwater ‘revap’ coefficient | 0.02–0.2 | Groundwater (.gw) | Replace |

GW_DELAY | Groundwater delay time (days) | 0–500 | Groundwater (.gw) | Replace |

^{a}SWAT parameters are constructed based on SWACUP user manual (Abbaspour 2013).

^{b}Prior distributions of aggregate parameters are based on Neitsch *et al.* (2011).

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.

### Model performance evaluation

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^{2} multiplied by the coefficient of the regression line (BR^{2}), and sum of the squares of the differences of the measured and simulated values after ranking (SSQR).

*et al.*2009); hence the modeled value represents variability and better fit to high and low flow events. where is the linear correlation coefficient between the observed flow and the simulated flow. is a measure of variability in the streamflow values (the standard deviation of the simulated flow over the standard deviation of the observed flow), and denotes the mean of simulated flow over the mean of observed flow (see Gupta

*et al.*2009).

^{2}(the coefficient of determination) multiplied by the coefficient of the regression line (Equation (10)). This function accounts for discrepancies in the magnitudes of two signals (depicted by

*b*) as well as their dynamics (depicted by R

^{2}): SSQR is similar to the mean square error method, and it is a measure to capture the fit of the frequency distribution to the observed and the simulated series. SSQR value of 0 indicates a perfect fit. The SSQ (Equation (11)) is defined as:

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.

## RESULTS AND DISCUSSION

### Model application and sensitivity analysis

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.

P–T | HG | P–M | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

SWAT parameter^{a} | P-value | t-stat | Rank | Optimal value | P-value | t-stat | Rank | Optimal value | P-value | t-stat | Rank | Optimal value |

v__CH_N2 | 0 | 7.48 | 1 | 0.10 | 0 | 4.19 | 2 | 0.101 | 0.144 | −1.46 | 7 | 0.10 |

r__SOL_K | 0.07 | 1.83 | 2 | 0.058 | 0.31 | 1 | 5 | 0.31 | 0.51 | −0.65 | 11 | −0.88 |

v__CH_K2 | 0.08 | 1.73 | 3 | 75.94 | 0.43 | −0.78 | 9 | 93.67 | 0.1 | 1.62 | 4 | 185.32 |

v__SLSUBBSN | 0.17 | 1.37 | 4 | 71.32 | 0.138 | 1.48 | 4 | 74.31 | 0 | 2.65 | 3 | 77.54 |

v__RCHRG_DP | 0.21 | −1.26 | 5 | 0.45 | 0.33 | −0.96 | 6 | 0.60 | 0.17 | 1.35 | 9 | 0.49 |

v__ESCO | 0.25 | −1.16 | 6 | 0.92 | 0 | 3.82 | 3 | 0.99 | 0 | 5.44 | 1 | 0.97 |

v__OV_N | 0.32 | −1 | 7 | 0.56 | 0.41 | −0.82 | 8 | 0.53 | 0.81 | −0.23 | 14 | 0.26 |

r__SOL_BD | 0.33 | 0.97 | 8 | 0.002 | 0.94 | 0.063 | 15 | 0.001 | 0.15 | 1.43 | 8 | 0.60 |

r__GWHT | 0.41 | 0.83 | 9 | 1.33 | 0.46 | 0.73 | 10 | 0.66 | 0.96 | −0.044 | 15 | −0.50 |

v__GW_REVAP | 0.44 | 0.78 | 10 | 0.087 | 0.48 | −0.7 | 11 | 0.18 | 0.142 | −1.46 | 6 | 0.08 |

v__EPCO | 0.5 | −0.63 | 11 | 0.55 | −0.17 | 0.86 | 14 | 0.17 | 0.706 | 0.37 | 13 | 0.92 |

r__SHALLST | 0.67 | 0.43 | 12 | −0.22 | 0.81 | −0.22 | 13 | 0.97 | 0.47 | 0.72 | 10 | −0.76 |

v__ALPHA_BF | −0.32 | 0.75 | 13 | 0.75 | 0.36 | −0.9 | 7 | 0.6 | 0.13 | −1.51 | 5 | 0.70 |

r__CN | −0.32 | 0.75 | 14 | 0.021 | 0 | 6.57 | 1 | 0.045 | 0 | 4.86 | 2 | 0.034 |

v__GW_DELAY | 0.94 | 0.07 | 15 | 207.71 | 0.7 | 0.37 | 12 | 307.99 | 0.69 | 0.4 | 12 | 249.86 |

P–T | HG | P–M | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

SWAT parameter^{a} | P-value | t-stat | Rank | Optimal value | P-value | t-stat | Rank | Optimal value | P-value | t-stat | Rank | Optimal value |

v__CH_N2 | 0 | 7.48 | 1 | 0.10 | 0 | 4.19 | 2 | 0.101 | 0.144 | −1.46 | 7 | 0.10 |

r__SOL_K | 0.07 | 1.83 | 2 | 0.058 | 0.31 | 1 | 5 | 0.31 | 0.51 | −0.65 | 11 | −0.88 |

v__CH_K2 | 0.08 | 1.73 | 3 | 75.94 | 0.43 | −0.78 | 9 | 93.67 | 0.1 | 1.62 | 4 | 185.32 |

v__SLSUBBSN | 0.17 | 1.37 | 4 | 71.32 | 0.138 | 1.48 | 4 | 74.31 | 0 | 2.65 | 3 | 77.54 |

v__RCHRG_DP | 0.21 | −1.26 | 5 | 0.45 | 0.33 | −0.96 | 6 | 0.60 | 0.17 | 1.35 | 9 | 0.49 |

v__ESCO | 0.25 | −1.16 | 6 | 0.92 | 0 | 3.82 | 3 | 0.99 | 0 | 5.44 | 1 | 0.97 |

v__OV_N | 0.32 | −1 | 7 | 0.56 | 0.41 | −0.82 | 8 | 0.53 | 0.81 | −0.23 | 14 | 0.26 |

r__SOL_BD | 0.33 | 0.97 | 8 | 0.002 | 0.94 | 0.063 | 15 | 0.001 | 0.15 | 1.43 | 8 | 0.60 |

r__GWHT | 0.41 | 0.83 | 9 | 1.33 | 0.46 | 0.73 | 10 | 0.66 | 0.96 | −0.044 | 15 | −0.50 |

v__GW_REVAP | 0.44 | 0.78 | 10 | 0.087 | 0.48 | −0.7 | 11 | 0.18 | 0.142 | −1.46 | 6 | 0.08 |

v__EPCO | 0.5 | −0.63 | 11 | 0.55 | −0.17 | 0.86 | 14 | 0.17 | 0.706 | 0.37 | 13 | 0.92 |

r__SHALLST | 0.67 | 0.43 | 12 | −0.22 | 0.81 | −0.22 | 13 | 0.97 | 0.47 | 0.72 | 10 | −0.76 |

v__ALPHA_BF | −0.32 | 0.75 | 13 | 0.75 | 0.36 | −0.9 | 7 | 0.6 | 0.13 | −1.51 | 5 | 0.70 |

r__CN | −0.32 | 0.75 | 14 | 0.021 | 0 | 6.57 | 1 | 0.045 | 0 | 4.86 | 2 | 0.034 |

v__GW_DELAY | 0.94 | 0.07 | 15 | 207.71 | 0.7 | 0.37 | 12 | 307.99 | 0.69 | 0.4 | 12 | 249.86 |

^{a}Aggregate parameters are constructed based on SWACUP user manual (Abbaspour 2013).

‘v__’ and ‘r__’ means a replacement, and a relative change to the initial parameter values, respectively.

### Effect of PET estimation models on streamflow simulation

*et al.*’s (2007) qualitative rank, the SWAT(P–T) model showed very good performance while both SWAT(HG) and SWAT(P–M) showed satisfactory performances in simulating the observed daily streamflow during 2003–2007. This result is in agreement with the findings of Lu

*et al.*(2005), Douglas

*et al.*(2009), and Rao

*et al.*(2011).

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.

*et al.*(2003, 2005) and Rao

*et al.*(2011) for the coastal plain watershed.

Quantiles | Observed (cm/s) | P–T (cm/s) | HG (cm/s) | P–M (cm/s) |
---|---|---|---|---|

Q10% | 1.11 | 10.83 | 16.08 | 10.04 |

Q25% | 8.01 | 14.92 | 21.09 | 18.72 |

Q50% | 23.13 | 23.58 | 28.42 | 27.82 |

Q75% | 46.16 | 47.46 | 43.65 | 44.65 |

Q90% | 85.40 | 84.32 | 79.78 | 75.26 |

Q95% | 124 | 110.58 | 105.06 | 96.51 |

Q99% | 173.70 | 163.90 | 175.98 | 149.86 |

Quantiles | Observed (cm/s) | P–T (cm/s) | HG (cm/s) | P–M (cm/s) |
---|---|---|---|---|

Q10% | 1.11 | 10.83 | 16.08 | 10.04 |

Q25% | 8.01 | 14.92 | 21.09 | 18.72 |

Q50% | 23.13 | 23.58 | 28.42 | 27.82 |

Q75% | 46.16 | 47.46 | 43.65 | 44.65 |

Q90% | 85.40 | 84.32 | 79.78 | 75.26 |

Q95% | 124 | 110.58 | 105.06 | 96.51 |

Q99% | 173.70 | 163.90 | 175.98 | 149.86 |

### Effect of PET estimation method on AET estimation

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.

### Spatial distribution of soil and groundwater physical properties

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.

## CONCLUSIONS

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.

## ACKNOWLEDGEMENTS

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.