Modeling the Fate and Transport of Malathion in the Pagsanjan-Lumban Basin , Philippines

Exposure to highly toxic pesticides could potentially cause cancer and disrupt the development of vital systems. Monitoring activities were performed to assess the level of contamination; however, these were costly, laborious, and short-term leading to insufficient monitoring data. However, the performance of the existing Soil and Water Assessment Tool (SWAT model) can be restricted by its two-phase partitioning approach, which is inadequate when it comes to simulating pesticides with limited dataset. This study developed a modified SWAT pesticide model to address these challenges. The modified model considered the three-phase partitioning model that classifies the pesticide into three forms: dissolved, particle-bound, and dissolved organic carbon (DOC)-associated pesticide. The addition of DOC-associated pesticide particles increases the scope of the pesticide model by also considering the adherence of pesticides to the organic carbon in the soil. The modified SWAT and original SWAT pesticide model was applied to the Pagsanjan-Lumban (PL) basin, a highly agricultural region. Malathion was chosen as the target pesticide since it is commonly used in the basin. The pesticide models simulated the fate and transport of malathion in the PL basin and showed the temporal pattern of selected subbasins. The sensitivity analyses revealed that application efficiency and settling velocity were the most sensitive parameters for the original and modified SWAT model, respectively. Degradation of particulate-phase malathion were also significant to both models. The rate of determination (R2) and Nash-Sutcliffe efficiency (NSE) values showed that the modified model (R2 = 0.52; NSE = 0.36) gave a slightly better performance compared to the original (R2 = 0.39; NSE = 0.18). Results from this study will be able to aid the government and private agriculture sectors to have an in-depth understanding in managing pesticide usage in agricultural watersheds.


Introduction
Agriculture has been substantial to the Philippine economy and has contributed 10.2% to 13.2% of the country's GDP in the past decade [1].To keep up with this demand, various kinds of pesticides were applied to different crops and vegetables to help increase food supplies and provide greater revenue for farmers.However, exposure to pesticides could potentially cause cancer and disrupt the development of vital systems (endocrine, reproductive, and immune systems) [2][3][4].Pesticide contamination in soil and water also has negative effects on the diversity of the flora and fauna of local areas thereby disturbing the existing ecosystem [5][6][7].Various kinds of pesticides are used in agriculture depending on the target pests.Hence, many kinds of chemicals exist and find their way into the groundwater, surface water, soils, and eventually drinking water [8][9][10][11].Several attempts have been made to monitor and map out their potential areas of contamination in the Philippines, especially in highly agricultural areas [12][13][14].
Laguna de Bay is the second largest freshwater lake in Southeast Asia and the largest in the Philippines.It is located east of Metro Manila, the Philippine capital, and is part of the Laguna de Bay basin.The basin has one of the fastest economic growth among others and it is a major water resource for agriculture, fisheries, and domestic use of the surrounding communities that has an estimated population of six million people [15].In the recent years, the lake has been threatened by waste discharges of the industrial, urban, and residential areas from the west and by intensive agricultural activities from the east [16].The presence of pesticides and other micropollutants led to the increasing levels of toxicity and fish-kill occurrences in the lake [17,18].Many efforts have been made to improve the water quality of the lake such as rehabilitation programs and cleanup operations within the vicinity of the basin.Monitoring activities were also performed to assess the level of contamination.However, these were short-term and limited to a few selections of pesticides [17][18][19].Applying environmental models to available monitoring datasets of pesticides can help broaden the understanding of the behavior of these micropollutants in the environment.However, existing modeling studies of Laguna de Bay basin lack watershed-scale analyses of pesticides used in its agricultural activities [14].
Processes driving pesticide fate and transport are on the whole well-known and are incorporated in various pesticide models operating at plot or watershed spatial scales, such as the crop model STICS (Simulateur mulTIdiscplinaire pour les Cultures Standard), MACRO (Water and solute transport in macroporous soils), PEARL (Pesticide Emission Assessment at Regional and Local scales), PRZM (Pesticide Root Zone Model), and SWAT (Soil and Water Assessment Tool) [20,21].Several studies have already shown that the watershed-scale SWAT model was an efficient tool to model pesticides fate and transport [22][23][24][25].However, the SWAT model performance can be restricted by its two-phase partitioning approach, which is inadequate when it comes to simulating pesticides with limited dataset.In this study, we modified the SWAT model by incorporating the three-phase partitioning model to improve the pesticide simulations, especially for watersheds with scarce dataset that are often common in developing countries.The modified model considered the three-phase partitioning model that classifies the pesticide into three forms: dissolved, particle-bound, and dissolved organic carbon (DOC)-associated pesticide.This approach is a first for pesticides and it differs from the original SWAT model that classified pesticides into two categories: pesticide sorbed into solid phase and pesticide in solution.The addition of DOC-associated pesticide particles increases the scope of the pesticide model by also considering the adherence of pesticides to the organic carbon in the soil.
We aimed to: (1) conduct a watershed-scale analysis of the fate and transport of pesticides, specifically malathion, and increase the accuracy of the simulated malathion loading using the modified pesticide model; and (2) perform a case study by applying the original SWAT model and modified pesticide model to a catchment with limited dataset, such as PL basin, and compare their performance.Malathion is an organophosphate insecticide used in PL basin for crops and vegetables.It is preferred by farmers due to its effectiveness against a wide range of pests and short half-life.The SWAT model was used to construct the watershed model for one of the subbasins of the Laguna de Bay basin, namely the Pagsanjan-Lumban (PL) basin.SWAT is a widely-used, physically-based hydrologic model that can predict the impact of water management practices [26,27].It can simulate the flowrate and the transport of nutrients, pesticides, and sediments in watersheds.Implementing a watershed-scale analysis of the fate and transport of the malathion in the Laguna de Bay basin using watershed models will give an insight on the dominant processes affecting pesticide loadings to the soil and water.Results from this study will aid the government and private agriculture sectors to have an in-depth understanding in managing pesticide usage in an agricultural watershed.[28].On the other hand, the western part has a Type III climate that has a short dry season, varying from 1 to 3 months, in December to February [28].Areas close to Mt. Banahaw at the southernmost part of the basin have relatively uniform rainfall distribution throughout the year [29].However, the basin in general experiences a dry period from November to April due to the rain shadow effect of the Sierra Madre mountain range while the wet period occurs for the remaining months [29].The average annual rainfall of the basin is 2996 mm, which mostly fall during the monsoon period.

Study Area
Figure 1  affecting pesticide loadings to the soil and water.Results from this study will aid the government and private agriculture sectors to have an in-depth understanding in managing pesticide usage in an agricultural watershed.

Study Area
PL basin is located at the southeastern part of Laguna de Bay basin in the Southern Tagalog Region (CALABARZON) of the Philippines.It has a catchment area of 454.45 km 2 (121°24′ E~121°37′ E, 14°37′ N~14°21′ N) that drains to Laguna de Bay.The watershed experiences two types of Philippine climate: (1) Type II; and (2) Type III.The eastern part of the basin experiences Type II climate that has no dry season with a very pronounced maximum rain period from December to February and a minimum rainfall period from March to May [28].On the other hand, the western part has a Type III climate that has a short dry season, varying from 1 to 3 months, in December to February [28].Areas close to Mt. Banahaw at the southernmost part of the basin have relatively uniform rainfall distribution throughout the year [29].However, the basin in general experiences a dry period from November to April due to the rain shadow effect of the Sierra Madre mountain range while the wet period occurs for the remaining months [29].The average annual rainfall of the basin is 2996 mm, which mostly fall during the monsoon period.
Figure 1    The basin has two major tributaries that branch out after the Lumban Station (outlet) as shown in Figure 1.Due to the presence of two reservoirs in the northern half of the basin, most of the discharge during dry days comes from the Pagsanjan River situated at the southern half of the PL basin.Pagsanjan River has a length of 54.1 km and a drainage area of 311.8 out of the 454.45 km 2 of PL basin, with a mean annual runoff of 53.1 m 3 •s −1 [29].The runoff pathways in the north were modified to collect the water in the reservoirs; hence, water only flows to the outlet during extreme rain events.

Monitoring Data
The locations of the monitoring stations for the flowrate, weather, and malathion are shown in Figure 1.[30,31], and the malathion concentrations were monitored at Lucban Station by Varca [16].A total of 26 sampling events at Lucban Station, with a frequency of at least two water samples a month, were carried out from December 2007 to November 2008 to measure the malathion concentrations.Each water sample was analyzed to measure the total concentration of malathion, which was used as comparison for the pesticide simulations in this study.

Hydrology Model
SWAT is a physically-based watershed model developed for the USDA Agricultural Research Service (ARS) to simulate the impact of land management practices on water, sediment, nutrients, and pesticide yields in large complex watersheds [32,33].The model operates at a daily time step and uses readily available inputs such as [34]: topography (DEM with a 90 m resolution from United States Geological Survey (USGS)/National Aeronautics and Space Administration Shuttle Radar Topography Mission), hydrography, weather data (INWARD and CFSR), landuse/land cover (USGS Global Land Cover Characterization database), and soil type (Food and Agriculture Organization).The delineation threshold of the PL basin was 2.5 km 2 , thus; it was divided into eight subbasins with 54 hydrological response units (HRU).Each of these HRUs is a unique combination of soil type, landuse, and slope.The threshold for soil, landuse, and slope was set to 0 to include non-agricultural areas in the basin that are less than 1% of the subbasin areas.
Table 1 summarizes the calibration and validation periods of the flowrate simulation.The calibration period was from September 2014 to May 2015 while the validation was from May 2014 to July 2014 and June 2016 to September 2016.The available flowrate dataset started from May 2014 until September 2016.We first compared the observed flowrate to the precipitation and noticed that the peaks of the flowrate did not match the precipitation for a period.This period was removed after concluding that the sensor was faulty at that time.The SWAT-Calibration and Uncertainty Program (SWAT-CUP) was then used to calibrate and validate the SWAT flowrate parameters shown in Table 2. Simulated flowrates in 2007 and 2008 were then applied to simulate malathion fate based on the available malathion dataset (December 2007 to November 2008).

Pesticide Modeling
The malathion loadings were then calibrated using the original SWAT and modified SWAT pesticide models.Both models need the management operation schedule to simulate the application of pesticide in the basin.In this study, malathion was applied to HRUs with "Tomato" (TOMA) as land cover.This is based on a previous study that summarized the pesticide usage of the farmers in the PL basin, which affects two HRUs from Subbasins 7 and 8 (shown in Figure 1) [35].TOMA was assumed as a collective representative of the vegetable crops in the PL basin that used malathion during the pesticide application period.Two planting seasons were implemented for TOMA.The first season starts in January while the second season starts in June.Malathion was then applied for 12 times for 3 months after planting at a rate of 0.57 kg/ha for every application (first season: January to March; second season: June to August).Based on this schedule, the pesticide models were then run to simulate malathion loading in the HRUs with TOMA as land cover.The models are further discussed in the next subsections.

Original SWAT Pesticide Model
The pesticide module in the SWAT model was applied to calibrate the malathion loadings in the PL basin.Figure 2 shows that the SWAT pesticide model used the two-phase partitioning approach that classify the pesticides as: pesticide sorbed into solid phase and pesticide in solution or liquid phase.Table 3 shows the pesticide parameters that describe the reaction and transport processes of malathion starting from the application (foliar, soil surface, and subsurface).These processes include degradation, infiltration, leaching, surface runoff, volatilization, and wash off mechanisms.

Pesticide Modeling
The malathion loadings were then calibrated using the original SWAT and modified SWAT pesticide models.Both models need the management operation schedule to simulate the application of pesticide in the basin.In this study, malathion was applied to HRUs with "Tomato" (TOMA) as land cover.This is based on a previous study that summarized the pesticide usage of the farmers in the PL basin, which affects two HRUs from Subbasins 7 and 8 (shown in Figure 1) [35].TOMA was assumed as a collective representative of the vegetable crops in the PL basin that used malathion during the pesticide application period.Two planting seasons were implemented for TOMA.The first season starts in January while the second season starts in June.Malathion was then applied for 12 times for 3 months after planting at a rate of 0.57 kg/ha for every application (first season: January to March; second season: June to August).Based on this schedule, the pesticide models were then run to simulate malathion loading in the HRUs with TOMA as land cover.The models are further discussed in the next subsections.

Original SWAT Pesticide Model
The pesticide module in the SWAT model was applied to calibrate the malathion loadings in the PL basin.Figure 2 shows that the SWAT pesticide model used the two-phase partitioning approach that classify the pesticides as: pesticide sorbed into solid phase and pesticide in solution or liquid phase.Table 3 shows the pesticide parameters that describe the reaction and transport processes of malathion starting from the application (foliar, soil surface, and subsurface).These processes include degradation, infiltration, leaching, surface runoff, volatilization, and wash off mechanisms.    Figure 3 shows the schematic diagram of the fate and transport of the pesticides for the modified SWAT pesticide model.The modified model applied in this study was based on the watershed-scale model from a previous study of the same authors about modeling the fate and transport of polycyclic aromatic hydrocarbons (PAH) and linking the PAH model with SWAT [36].This study further developed the model to include the pesticide application based on the original SWAT model and other equations related to the fate and transport of pesticides.The approach of the original 2-phase partitioning SWAT model on the fate and transport of pesticides was modified by considering the three-phase partitioning model shown in Figure 2. Figure 4 shows the diagram of using MATLAB as platform for the modified approach.Table 4 shows the parameters of the modified pesticide model.Figure 3 shows the schematic diagram of the fate and transport of the pesticides for the modified SWAT pesticide model.The modified model applied in this study was based on the watershed-scale model from a previous study of the same authors about modeling the fate and transport of polycyclic aromatic hydrocarbons (PAH) and linking the PAH model with SWAT [36].This study further developed the model to include the pesticide application based on the original SWAT model and other equations related to the fate and transport of pesticides.The approach of the original 2-phase partitioning SWAT model on the fate and transport of pesticides was modified by considering the three-phase partitioning model shown in Figure 2. Figure 4 shows the diagram of using MATLAB as platform for the modified approach.Table 4 shows the parameters of the modified pesticide model.This study further developed the model to include the pesticide application based on the original SWAT model and other equations related to the fate and transport of pesticides.The approach of the original 2-phase partitioning SWAT model on the fate and transport of pesticides was modified by considering the three-phase partitioning model shown in Figure 2. Figure 4 shows the diagram of using MATLAB as platform for the modified approach.Table 4 shows the parameters of the modified pesticide model.Temperature adjustment factor for DOC-associated pesticide -

Accumulation of Pesticide
Pesticides are usually distributed to the foliage and soil surface during application.In this case, separate equations were applied to calculate the amount of pesticides on the foliage and the soil surface [37].This approach is similar to the original SWAT model.The amount of pesticides on foliage (pst f , kg of pesticide ha −1 ) and on the soil surface (pst surf , kg of pesticide ha −1 ) were determined by the equations below [37]: where gc is the fraction of the ground surface covered by plants (-) and pst is the efficient amount of pesticide applied (kg of pesticide ha −1 ).These terms were determined by the following equations [37]: where erfc is the complementary error function (-), LAI is the leaf area index (-), ap ef is the pesticide application efficiency (-), and pst is the original amount of pesticide applied (kg of pesticide ha −1 ).Pesticide on the foliage was assumed to be affected by wash off due to rain events and degradation.The amount of pesticides washed off by precipitation from the plants (pst f,wsh , kg of pesticide ha −1 ) were determined using the equation below [37]: where fr wsh is the wash off fraction for the pesticide on the foliage (-).Degradation of pesticide on the foliage was then calculated using the following equation [37]: where pst f,t is the amount of pesticide on the foliage at time t (kg of pesticide ha −1 ), pst f,0 is the initial amount of pesticide on the foliage (kg of pesticide ha −1 ), k p,f is the degradation rate constant of the pesticide (day −1 ), and t is time (day).

Pesticides on Soil
The original SWAT model assumed that the pesticide is either sorbed to the solid phase or dissolved in solution [37].The three-phase partitioning model was used to estimate the different forms of pesticide in the soil.This method classifies the pesticide into three classes: dissolved pesticide, pesticide adsorbed on dissolved organic carbon ([DOC]), and pesticide adsorbed on particles.The total pesticide and dissolved pesticide concentrations in the bulk saturated soil are defined by the following equations [38]: where C bs t is the total pesticide concentration in the bulk soil (kg•L −1 bulk soil), C bs p is the concentration of particle-bound pesticide (kg•L −1 bulk soil), C bs d is the dissolved pesticide concentrations in the bulk saturated soil (kg•L −1 bulk soil), C bs fd is the dissolved pesticide, and C bs DOC is the DOC-associated pesticide.Combining Equations ( 7) and ( 8) yields the following equation: The terms C bs p and C bs D•C can also be determined by multiplying C bs fd by the coefficients shown in Equations ( 10)-( 12) below [38]: [DOC] = f DOC × OM where r sw is the soil-to-water ratio (kg•L −1 ), K sw is the soil-water distribution coefficient [L•kg −1 ], K DOC is the dissolved organic carbon-water partition coefficient (L•kg −1 ), f DOC is the fraction of DOC, and OM is the concentration of the organic matter (kg•L −1 ).OM was calculated by dividing the mass of the soil carbon in the soil organic matter with the water yield, which were both simulated by the original SWAT model [36].C bs p , C bs fd , and C bs D•C from Equation ( 9) can also be determined by using the pesticide in bulk soil fractions [36,38]: where f p (-), W cp (kg•L −1 bulk soil), and µ p (s −1 ) are the fraction, wash off load, and rate constant of particle-bound pesticide, ε pstsed is the enrichment ratio (-), f d fd (-), Wcf (kg•L −1 bulk soil), and µ f (s −1 ) are the fraction, wash off load, and rate constant of dissolved pesticide, and f d DOC (-), W cf (kg•L −1 bulk soil), and µ DOC (s −1 ) are the fraction, wash off load, and rate constant of DOC-bound pesticide.The enrichment ratio was calculated using this equation [37] : where x 1 and x 2 are the enrichment ratio coefficients, sed is the sediment yield (metric tons), and WY is the water yield (mm H 2 O).Surface runoff of the suspended particles to the channel was also considered in the wash off loads and wash off fractions of the pesticide.The suspended particles affected by the wash off mechanism are described in the following equations [36]: where W cp , W cf , and W cDOC are the wash off loads of the particle-bound, dissolved, and DOC-associated PAHs exported to the river via runoff (kg•L −1 bulk soil); C p1 , C f1 , and C DOC1 are the wash off coefficients for the particle-bound, dissolved, and DOC-associated PAHs (-); q is the runoff rate per unit area; and C p2 , C f2 , and C DOC2 are the wash off exponents for the particle-bound, dissolved, and DOC-associated PAHs (-), respectively.The degradation term, previously shown in Equations ( 13)-( 15), is a temperature-dependent rate constant expressed as the first-order reaction.The rate constant for the three phases was defined as: where µ i is the initial rate constant for the pesticide [s −1 ], θ is the temperature adjustment factor for pesticide (-), and T is the temperature ( • C). µ and θ were then calibrated for the particle-bound, dissolved, and DOC-associated pesticide.Pesticide loadings on the soil and into the water were computed by [38]: where C p and C w are the pesticide loading on soil [pesticide per solid mass] and in water [pesticide per fluid mass], respectively, and C wfinal is the final pesticide loading in water.

Pesticides in Water
The pesticides in the channels are then subjected to the following mechanisms and processes due to water movement and reactivity of the pesticide with other components present in the water: advection, dispersion, photodegradation, and settling processes upon entering the channel [36].The concentration of the pesticide in the waterbody was determined after considering these processes.This is expressed by the modified advection-dispersion equation: where C is the concentration of pesticide in water (g•m −3 ), t is time (day), x is distance (m), u is the velocity of the water (m•s −1 ), D L is the dispersion coefficient (m 2 •d −1 ), f is the fraction of the particulate pesticide in water (-), v s is the settling velocity (m•s −1 ), h is the depth of the channel (m), a is the photodegradation coefficient (m 2 •MJ −1 •d −1 ), and I is the solar intensity (MJ•m −2 ).The photodegradation term in the equation was added to the advection-dispersion equation to fit the pesticide model.

Sensitivity Analyses
The sensitivity analyses for the flowrate and pesticide were simultaneously done with the calibration.SWAT-CUP applies the Latin-Hypercube (LH) sampling method, which is based on the Monte Carlo simulation, and set the Nash-Sutcliffe efficiency (NSE) as the objective function (Section 2.6).This is a robust method that requires a large number of simulations and computational resources [39].LH sampling randomly assigns a value within the permitted range of each parameter to complete a parameter set.The One-factor-At-a-Time (OAT) sensitivity test was then commenced after the sampling.This method takes one parameter for each run and changes its value to determine how each parameter affects the results.Sensitive parameters were determined based on their probability values or p-values.Parameters with less than 0.01 p-value were labeled as sensitive.The LH-OAT method was applied to the SWAT hydrology and pesticide models in this study.

Evaluation Criteria
The coefficient of determination (R 2 ), Nash-Sutcliffe Efficiency coefficient (NSE), root mean square error (RMSE), and percent bias (PBIAS) were then calculated to evaluate the SWAT and the modified model.These statistical indices can determine whether the model performance is satisfactory or unsatisfactory [40].R 2 and NSE were calculated for the calibration and validation periods of the flowrate (Table 1), while R 2 , NSE, RMSE and PBIAS were determined for pesticide.For a daily time-step, the model is acceptable or satisfactory when the R 2 and NSE values are greater than 0.5 and when the RMSE and PBIAS values reach the optimal value of 0. A lower RMSE value is an indicator that the model has less residual variance [41].PBIAS shows if the model has under-or overestimated the results compared to the observations, and lower absolute PBIAS values indicate more accurate model simulations [42].
SWAT-CUP also has other criteria to quantify the strength and uncertainties of the calibration analysis.Given the small dataset of this study, the P-factor and R-factor of the iteration were also noted.P-factor represents percentage of the simulated data covered by the 95% uncertainty band (95PPU) and it ranges from 0 to 1 [43].A P-factor value greater than 0.5 is desirable since it entails that more than 50% of the simulated data are acceptable.R-factor, on the other hand, estimates strength of the calibration by dividing the average thickness of the 95PPU band with the standard deviation of the measured data [43].The range of the R-factor starts from 0 to infinity; high and low R-factor values correspond to a thicker and thinner 95 PPU bands, respectively.A P-factor of 1.0 and an R-factor of 0 signifies a perfect simulation of the observed data.

Flowrate Calibration
Figure 5 shows the comparison of the observed and SWAT-simulated flowrate at Lumban Station.The model did not implement the presence of reservoirs and paddy fields due to the lack of information; hence, it was assumed that these factors would affect the flowrate calibration when it comes to water storage.The calibration process yielded an R 2 value of 0.42 and an NSE of 0.22 while the validation period has 0.10 and −2.87, respectively.These values are less than the desired value of 0.5 [40], showing that the SWAT model is slightly unfit and underperforming.The lack of information regarding the management of the two reservoirs (storage, release, and distribution) may have significantly degraded the performance of the model [44].The SWAT model of the PL basin was limited to a small dataset, which was a major disadvantage during the calibration (227 observed flowrate) and validation (158 observed flowrate) of the flowrate.The intense rainfall events happened in the calibration period; hence, the entire basin drains to the outlet, as mentioned in Section 2.1.This is in contrast with the validation period that included the days with dry to mild rainfall, which only drains the southern part of the basin.The 95PPU band of the flowrate calibration in Figure 5 has a P-factor value greater than 0.5 and an R-factor value of less than 1.This result indicates that more than 50% of the simulated flowrate is within the acceptable uncertainty [43].The validation period, on the other hand, has a P-factor of 0.18 and an R-factor of 0.72, indicating that only 18% of the simulated flowrate during this period is within the uncertainty range.
Nine parameters were found to have a significant effect on the flowrate (Table 5).Surface runoff lag coefficient (SURLAG) was the most sensitive parameter with a calibration value of 0.051 h.This value is negligible compared to the default value in SWAT, which is four hours.As SURLAG decreases, more water is stored in the basin, indicating that surface runoff does not go directly to the channel in the PL basin and is instead stored elsewhere [45].The presence of two reservoirs in the northern half of the basin, mentioned in Section 2.1, may have affected the SURLAG value.A huge percentage of the runoff was immediately collected by the reservoirs instead of getting released across the basin.The water can also be stored in the foliage and paddies that are present in the basin.This is supported by the other sensitive parameters: Manning's "n" value for the tributary channels (CH_N1), second most sensitive; the initial SCS runoff curve number for moisture condition II (CN2), third; baseflow alpha factor (ALPHA_BF), fourth; and effective hydraulic conductivity in tributary channel alluvium (CH_K1), fifth.CN2 signifies the soil permeability of the basin [46].Increasing CN2 values is associated with the increase of imperviousness of the basin.This can be related to how urbanized the basin is.During calibration, SWAT-CUP relatively changes the CN2 value at the HRU level since HRUs have varying CN2 values.In general, CN2 in PL basin has a calibrated value of −0.096 or −9.6%, indicating a general decrease in the CN2 values of each HRU.Hence, the basin is slightly more pervious compared to the default values suggested by SWAT.This result was expected since PL basin is highly agricultural [29].ALPHA_BF, on the other hand, is the baseflow recession constant, a direct index of groundwater flow response and is a basin-wide parameter.It has a calibrated value of 0.46, suggesting an intermediate or average response to the change in recharge, slightly leaning towards the slow response (slow response: 0.1-0.3;rapid response: 0.9-1.0).The result can indicate that it is possible for water to be stored in the shallow aquifer that can also be associated to the CH_K1 value.CH_K1 controls the transmission losses from the surface runoff in the tributary.It is determined by the type of bed materials present in the channel bed.The calibrated result of 48.4 mm•h −1 falls on the moderately high loss rate that ranges from 25 mm•h −1 to 76 mm•h −1 .This value characterizes the  CN2 signifies the soil permeability of the basin [46].Increasing CN2 values is associated with the increase of imperviousness of the basin.This can be related to how urbanized the basin is.During calibration, SWAT-CUP relatively changes the CN2 value at the HRU level since HRUs have varying CN2 values.In general, CN2 in PL basin has a calibrated value of −0.096 or −9.6%, indicating a general decrease in the CN2 values of each HRU.Hence, the basin is slightly more pervious compared to the default values suggested by SWAT.This result was expected since PL basin is highly agricultural [29].ALPHA_BF, on the other hand, is the baseflow recession constant, a direct index of groundwater flow response and is a basin-wide parameter.It has a calibrated value of 0.46, suggesting an intermediate or average response to the change in recharge, slightly leaning towards the slow response (slow response: 0.1-0.3;rapid response: 0.9-1.0).The result can indicate that it is possible for water to be stored in the shallow aquifer that can also be associated to the CH_K1 value.CH_K1 controls the transmission losses from the surface runoff in the tributary.It is determined by the type of bed materials present in the channel bed.The calibrated result of 48.4 mm•h −1 falls on the moderately high loss rate that ranges from 25 mm•h −1 to 76 mm•h −1 .This value characterizes the tributary bed material as sand and gravel mixture with low silt-clay content, which makes it easier for the transmission losses of to percolate into the shallow aquifer.Another sensitive parameter that has a similar description to CH_K1 is the effective hydraulic conductivity in the main channel alluvium (CH_K2).The calibrated value of CH_K2, 42.2 mm•h −1 , falls within the same range as CH_K1, indicating that the main channel has similar characteristics as the tributary channel.Having a lower CH_K2 value compared to CH_K1 signifies that the bed material of the main channel is a slightly more consolidated compared to the tributary channels.The second most sensitive parameter also describes the tributary channel in the basin, CH_N1.CH_N1 has a calibrated value of 0.12, which indicates the channel is well maintained and full of weeds and brushes (excavated or dredged) or it is heavy timbered with lots of vegetation as well (natural stream).The deep aquifer percolation fraction (RCHRG_DP) was also found to be sensitive, with a calibrated value of 0.97.This result indicates that a huge fraction of percolation from the root zone recharges the deep aquifer [47].However, the value is extremely high compared to the other studies with high RCHRG_DP values.Schuol et al. [48] estimated a range of 0.4 to 0.65 for the West Africa subcontinent, while Me et al. [47] yielded a value of 0.87 for a New Zealand catchment with mixed landuse.Though it can be assumed that the RCHRG_DP value for PL basin is also high, it should also be noted that this parameter may have been affected by the discrepancies formulated from the limited dataset (Section 2.2) and that the fraction (0.97) is too high for this basin.Lastly, the average slope length (SLSUBBSN) and Manning's value for overland flow (OV_N) have calibrated values of 75.8 m and 0.33, respectively.

Pesticide Calibration with SWAT Pesticide Model
After the flow calibration and validation, SWAT-CUP was applied to calibrate and analyze the sensitive parameters of the SWAT pesticide model.The calibration processes of the flowrate and pesticide were done separately to incorporate the same calibrated values of the flowrate parameters for the original SWAT model and the modified model.Figure 6 shows the observed and the simulated results of the malathion concentrations at Lucban Station using the original SWAT model.The calibration process yielded an R 2 of 0.39, an NSE of 0.18, a PBIAS value of 59.2%, and an RMSE of 3492.23 mg.Five parameters were found to be significant in the SWAT pesticide model (Table 6).The most sensitive parameter was the application efficiency (AP_EF) of the malathion with a calibrated value of 0.13.This indicates that only 13% of the applied malathion is deposited on the foliage and soil surface while the rest are lost in the atmosphere, which can be due to the type of pesticide application and management practices of the farmers in the basin.The pesticide partition coefficient between water and sediment in the reach (CHPST_KOC) was the second most sensitive parameter with a fitted value of 0.0012.A low value indicates that malathion is highly mobile in the water in its dissolved form, thereby increasing its potential for long-distance transport [16,49].This parameter is followed by the degradation half-life of malathion on the soil surface (HLIFE_S), the reaction coefficient in the channel (CHPST_REA), and the soil adsorption coefficient for soil organic carbon (SKOC) that have calibrated values of 6.26 day −1 , 0.037 day −1 , and 3594.52L•kg −1 , respectively.The sensitivity analysis determined that the pesticide application, degradation in the soil, and the sediment interaction and reactivity of malathion in water were the important processes that influenced the fate and transport of malathion in the basin.

Pesticides Calibration with Modified Pesticide Model
Table 7 revealed that the settling velocity (vs) was the most sensitive parameter, which suggest that the interaction of malathion particles in the water greatly affects the malathion transport.Degradation rate constant of the dissolved malathion particles (µk,fd) was the second most sensitive, which can be attributed to how malathion readily dissolves in water compared to soil.This was followed by, in no particular order, the wash off coefficient (Cfd,1) and exponent (Cfd,2), diffusion coefficient (En), porosity (poro), soil density (ρsoil), organic carbon fraction in soil (foc), and temperature adjustment factor of the particle-bound malathion (θk,p).Malathion parameters that are associated with the particulate phase were sensitive for both models.However, the specific parameters are not exactly the same.The modified model has a more detailed formalism since it applied the three-phase partitioning model.This gave a visual understanding of the different forms of malathion that were greatly affected by wash off, which is the dissolved malathion.In this case, it can be assumed that dissolved malathion is more susceptible to wash off and most likely to end up in the channel compared to the other two malathion phases, particle-bound and DOC-associated malathion.Aside from the malathion-specific parameters, soil properties were also important such as the soil density, porosity, and organic carbon fraction of the soil.

Comparison of Observed and Simulated Malathion Loading
Figure 7 shows the observed and simulated values of the monitoring dataset.The modified model has 0.52, 0.36, 48.6%, and 3088.05mg values for R 2 , NSE, PBIAS, and RMSE, respectively.Based on these evaluation criteria, the modified model performed better compared to the SWAT model (statistics mentioned in Section 3.2).

Pesticides Calibration with Modified Pesticide Model
Table 7 revealed that the settling velocity (v s ) was the most sensitive parameter, which suggest that the interaction of malathion particles in the water greatly affects the malathion transport.Degradation rate constant of the dissolved malathion particles (µ k,fd ) was the second most sensitive, which can be attributed to how malathion readily dissolves in water compared to soil.This was followed by, in no particular order, the wash off coefficient (C fd,1 ) and exponent (C fd,2 ), diffusion coefficient (En), porosity (poro), soil density (ρ soil ), organic carbon fraction in soil (f oc ), and temperature adjustment factor of the particle-bound malathion (θ k,p ).Malathion parameters that are associated with the particulate phase were sensitive for both models.However, the specific parameters are not exactly the same.The modified model has a more detailed formalism since it applied the three-phase partitioning model.This gave a visual understanding of the different forms of malathion that were greatly affected by wash off, which is the dissolved malathion.In this case, it can be assumed that dissolved malathion is more susceptible to wash off and most likely to end up in the channel compared to the other two malathion phases, particle-bound and DOC-associated malathion.Aside from the malathion-specific parameters, soil properties were also important such as the soil density, porosity, and organic carbon fraction of the soil.

Comparison of Observed and Simulated Malathion Loading
Figure 7 shows the observed and simulated values of the monitoring dataset.The modified model has 0.52, 0.36, 48.6%, and 3088.05mg values for R 2 , NSE, PBIAS, and RMSE, respectively.Based on these evaluation criteria, the modified model performed better compared to the SWAT model (statistics mentioned in Section 3.2).The observed data were compared to the simulated loading by SWAT and the modified model in Figure 8.The figure includes the logarithmic (Figure 8a) and actual (Figure 8b) scale of the values, which reveals the similarities and differences between the models.Both models were able to achieve a small deviation between the observed and simulated low malathion loading.Comparing the two simulations, the modified model captured more low values compared to SWAT.However, the models were poorly able to simulate the high values as seen in the two plots (Figure 7), the results were comparable for both models.Figure 7 shows the time series of the malathion loading simulated by the SWAT and modified models compared to the SWAT-projected flowrate from January 2007 to December 2008 at the Lucban Station.The malathion simulations peaked during the duration of the pesticide application.However, the peaks of modified model showed more consistency compared to the increasing peaks of the SWAT model.Both models have similar peak levels at the fourth peak, but the SWAT model gave a more distinct pattern compared to the modified model.The observed data were compared to the simulated loading by SWAT and the modified model in Figure 8.The figure includes the logarithmic (Figure 8a) and actual (Figure 8b) scale of the values, which reveals the similarities and differences between the models.Both models were able to achieve a small deviation between the observed and simulated low malathion loading.Comparing the two simulations, the modified model captured more low values compared to SWAT.However, the models were poorly able to simulate the high values as seen in the two plots (Figure 7), the results were comparable for both models.Figure 7 shows the time series of the malathion loading simulated by the SWAT and modified models compared to the SWAT-projected flowrate from January 2007 to December 2008 at the Lucban Station.The malathion simulations peaked during the duration of the pesticide application.However, the peaks of modified model showed more consistency compared to the increasing peaks of the SWAT model.Both models have similar peak levels at the fourth peak, but the SWAT model gave a more distinct pattern compared to the modified model.The observed data were compared to the simulated loading by SWAT and the modified model in Figure 8.The figure includes the logarithmic (Figure 8a) and actual (Figure 8b) scale of the values, which reveals the similarities and differences between the models.Both models were able to achieve a small deviation between the observed and simulated low malathion loading.Comparing the two simulations, the modified model captured more low values compared to SWAT.However, the models were poorly able to simulate the high values as seen in the two plots (Figure 7), the results were comparable for both models.Figure 7 shows the time series of the malathion loading simulated by the SWAT and modified models compared to the SWAT-projected flowrate from January 2007 to December 2008 at the Lucban Station.The malathion simulations peaked during the duration of the pesticide application.However, the peaks of modified model showed more consistency compared to the increasing peaks of the SWAT model.Both models have similar peak levels at the fourth peak, but the SWAT model gave a more distinct pattern compared to the modified model.

Conclusions
Insufficient data were one of the limitations of this study that imply high uncertainty in both models.However, building a SWAT model for the flowrate calibration was made possible in this study.Equations that are relevant to the fate and transport of pesticides were added to the modified model prior to simulations.The objectives of this study were also met by simulating the malathion loading in the PL basin using the SWAT and modified model.Hence, the following conclusions were derived: 1.
The sensitivity analysis of the hydrology model revealed that the flowrate of the PL basin is greatly influenced by the perviousness of the soil and the characteristics of the tributary channel that stores and retains the water in the basin.

2.
The modified pesticide model gave a slightly better performance compared to the original SWAT model, considering the statistical analyses performed (R 2 , NSE, PBIAS, and RMSE).

3.
Application efficiency was the most sensitive parameters for the original SWAT model, suggesting a possible need to improve the pesticide application and management practices of farmers in the basin, while settling velocity was the most sensitive for the modified models.Parameters associated with particulate-phase malathion, especially the degradation of particle-bound malathion, were also significant to both models.4.
The temporal patterns of the target subbasin simulated by the models showed that the modified model has more consistent peaks during the duration of pesticide application compared to SWAT.
This study focused on the comparison of the outcomes of the modified model with the commonly used SWAT hydrological model.The modified model and the original SWAT model were able to identify similar sensitive parameters.However, further development of the model is needed to incorporate pesticide application scenarios and interaction of soil and water media to the atmosphere.
PL basin is located at the southeastern part of Laguna de Bay basin in the Southern Tagalog Region (CALABARZON) of the Philippines.It has a catchment area of 454.45 km 2 (121 • 24 E~121 • 37 E, 14 • 37 N~14 • 21 N) that drains to Laguna de Bay.The watershed experiences two types of Philippine climate: (1) Type II; and (2) Type III.The eastern part of the basin experiences Type II climate that has no dry season with a very pronounced maximum rain period from December to February and a minimum rainfall period from March to May shows the Digital Elevation Model (DEM) of the PL basin.The areas near Mt.Banahaw at the southern region have the highest elevation, ranging from 560 m to 2170 m, while the eastern region near Sierra Madre ranges from 350 m to 560 m.The region near the outlet, including Lumban delta, has the lowest elevation, which ranges from 0 m to 200 m.Negative values can also be observed within the Lumban delta indicating that the elevations are below sea level and are often submerged in water.The outlet of the basin was set at the Lumban Station before the river branched out to the Lumban delta to exclude the possibility of water intrusion from the lake.Water 2017, 9, 451 3 of 18 shows the Digital Elevation Model (DEM) of the PL basin.The areas near Mt.Banahaw at the southern region have the highest elevation, ranging from 560 m to 2170 m, while the eastern region near Sierra Madre ranges from 350 m to 560 m.The region near the outlet, including Lumban delta, has the lowest elevation, which ranges from 0 m to 200 m.Negative values can also be observed within the Lumban delta indicating that the elevations are below sea level and are often submerged in water.The outlet of the basin was set at the Lumban Station before the river branched out to the Lumban delta to exclude the possibility of water intrusion from the lake.

Figure 1 .
Figure 1.The digital elevation model of the Pagsanjan-Lumban watershed.

Figure 1 .
Figure 1.The digital elevation model of the Pagsanjan-Lumban watershed.

Figure 2 .Table 3 .CHPST_KOC. 1 CHPST_REA 10 CHPST_STL
Figure 2.Original SWAT model applies the two-phase partitioning approach: pesticide in liquid phase and pesticide sorbed to the solid phase.Modified pesticide model assumes the three-phase partitioning model: dissolved pesticide, particle-bound pesticide, and DOC-associated pesticide.

Figure 2 .
Figure 2.Original SWAT model applies the two-phase partitioning approach: pesticide in liquid phase and pesticide sorbed to the solid phase.Modified pesticide model assumes the three-phase partitioning model: dissolved pesticide, particle-bound pesticide, and DOC-associated pesticide.

Figure 3 .
Figure 3. Fate and transport diagram of pesticides in the environment with the modified SWAT model.The three-phased partitioning model approach that we applied for the pesticide was also based on a previous study of the same authors about modeling the fate and transport of polycyclic aromatic hydrocarbons (PAH) and linking the PAH model with SWAT.

Figure 4 .
Figure 4. Flow diagram of the modified pesticide model.SWAT output and management plan of malathion were used as input in MATLAB for the pesticide simulations.

Figure 3 .Figure 3
Figure 3. Fate and transport diagram of pesticides in the environment with the modified SWAT model.The three-phased partitioning model approach that we applied for the pesticide was also based on a previous study of the same authors about modeling the fate and transport of polycyclic aromatic hydrocarbons (PAH) and linking the PAH model with SWAT.

Figure 3 .
Figure 3. Fate and transport diagram of pesticides in the environment with the modified SWAT model.The three-phased partitioning model approach that we applied for the pesticide was also based on a previous study of the same authors about modeling the fate and transport of polycyclic aromatic hydrocarbons (PAH) and linking the PAH model with SWAT.

Figure 4 .
Figure 4. Flow diagram of the modified pesticide model.SWAT output and management plan of malathion were used as input in MATLAB for the pesticide simulations.

Figure 4 .
Figure 4. Flow diagram of the modified pesticide model.SWAT output and management plan of malathion were used as input in MATLAB for the pesticide simulations.

Figure 5 .
Figure 5. Observed and simulated flowrate (m 3 •s −1 ) at Lumban Station over the calibration (September 2014-May 2015) and the validation (May 2014-July 2014 and June 2016-September 2016) periods, together with the 95% uncertainty band (95 PPU).The flowrate is plotted against the precipitation (mm) for the whole period.

Figure 5 .
Figure 5. Observed and simulated flowrate (m 3 •s −1 ) at Lumban Station over the calibration (September 2014-May 2015) and the validation (May 2014-July 2014 and June 2016-September 2016) periods, together with the 95% uncertainty band (95 PPU).The flowrate is plotted against the precipitation (mm) for the whole period.

Figure 7 .
Figure 7. (a) Logarithmic scale (-); and (b) normal scale (milligram) of the observed and simulated malathion loadings at Lucban Station.This shows the comparison of the observed malathion data and simulated loading by the original SWAT and modified SWAT model.

Figure 8 .
Figure 8.Time series of malathion loading simulated by the original SWAT and the modified SWAT models, observed malathion loading at the Lucban Station, and SWAT simulated flowrate.

Figure 7 .
Figure 7. (a) Logarithmic scale (-); and (b) normal scale (milligram) of the observed and simulated malathion loadings at Lucban Station.This shows the comparison of the observed malathion data and simulated loading by the original SWAT and modified SWAT model.

Water 2017, 9 , 451 15 of 18 Figure 7 .
Figure 7. (a) Logarithmic scale (-); and (b) normal scale (milligram) of the observed and simulated malathion loadings at Lucban Station.This shows the comparison of the observed malathion data and simulated loading by the original SWAT and modified SWAT model.

Figure 8 .
Figure 8.Time series of malathion loading simulated by the original SWAT and the modified SWAT models, observed malathion loading at the Lucban Station, and SWAT simulated flowrate.

Figure 8 .
Figure 8.Time series of malathion loading simulated by the original SWAT and the modified SWAT models, observed malathion loading at the Lucban Station, and SWAT simulated flowrate.
The daily flowrate data from May 2014 to October 2016 at Lumban Station and the weather data (precipitation, temperature, humidity, and wind speed) at Cavinti Station from 2014 to 2016 were acquired from the Integrated National Watershed Research and Development Project (INWARD), the weather data from 1979 to 2014, including the solar values, were generated from the Climate Forecast System Reanalysis (CFSR) of the Global Weather Data for SWAT

Table 1 .
Summary of the calibration and validation periods.

Table 2 .
Soil and Water Assessment Tool (SWAT) flowrate parameters for calibration and sensitivity analysis.

Table 3 .
Pesticide parameters for calibration and sensitivity analysis of the original SWAT model.

Table 4 .
Pesticide parameters for calibration and sensitivity analysis of the modified model.

Table 5 .
SWAT parameters for calibration and sensitivity analysis of the original SWAT pesticide model.

Table 5 .
SWAT parameters for calibration and sensitivity analysis of the original SWAT pesticide model.

Table 6 .
Sensitive parameters of the SWAT pesticide model.

Table 7 .
Sensitive parameters of the modified pesticide model.

Table 7 .
Sensitive parameters of the modified pesticide model.