Spatial Model Assessment of P Transport from Soils to Waterways in an Eastern Mediterranean Watershed

P index is a management tool commonly used to identify critical source areas (CSAs) in agro-catchments. We tested the applicability of several P-index models adjusted to Eastern Mediterranean conditions. On the basis of model structure and data requirements, we selected the Arkansas model and two models with the RUSLE equation and runoff curve number (RCN). Concurrently, we developed a GIS-based Hermon-P model which was designed to simulate rainfall–runoff events representing the major nutrient-transport mechanism in Eastern Mediterranean. The P index values computed by the Arkansas and RUSLE models exhibited low correlation (r 2 < 0.32) with the measured soluble reactive (SRP) and total P (TP), while the RCN model result correlations were somewhat higher (r 2 = 0.53 for SRP and 0.45 for TP). High correlations between the calculated and measured P during rainfall–runoff events were only achieved with the Hermon model (r 2 = 0.77 to 0.9). These high coefficients resulted from avoiding subjective categorization of the continuous variables and using the measured site-specific erosional predictors instead. On one occasion, during the first significant runoff event of the year, the Hermon model failed to predict total P in the stream water (r 2 = 0.14) because of considerable resuspension of stream sediments. Most of the P-index models are based on the perceptual transfer-continuum framework that was developed for temperate agro-catchments; this framework does not consider P resuspension along streams during 263 rainfall–runoff events. Hence, a new set of equations should be added to the P index to account for potential resuspension in Eastern Mediterranean streams at the beginning of the hydrological year.


Introduction
Nutrient transport from farmland to waterways may lead to eutrophication in downstream water bodies followed by massive algal growth and a significant decrease in water quality [1,2].Algal composition and bloom intensity are influenced by N and P concentrations, temperature, light, and water-column stability, among other physicochemical parameters [3,4].Since N-fixing cyanobacteria exhibit a clear advantage over non-fixing bacteria in lakes with low N:P ratio, P has become the limiting nutrient in many aqueous environs [3,5,6].Sharpley and Rekolainen [7] and Howarth et al. [8] demonstrated a clear linkage between water quality in lakes and natural processes and management practices in upstream watersheds.Several studies have stressed the importance of protecting the hydroecological system in watersheds to minimize the transport of nutrients, particularly P. Since P transport from the terrestrial system to waterways is strongly affected by spatiotemporal factors, it is important to determine the relative contribution of the major attributes that control P transport in a given watershed [9].However, Gburek et al. [10] showed that only a small portion of the watershed contributes a significant amount of P, which they termed critical source area (CSA).They found that most of the P in a small agro-watershed arrives from areas adjacent to the streams (~60 m).In support of this notion, McDowell [11] and James et al. [12] showed that cattle grazing in areas near the stream were the major contributing factor in P transport to the tributary, while little P is generated from grazing in more distant pastures.In the Mediterranean climate, P fate and transport is highly dependent on a few rainfall events that may generate significant runoff episodes [13,14].To minimize P transport in Mediterranean agro-watersheds, it is imperative to identify the CSAs and evaluate their temporal nature.To achieve this, a spatial model that will run over several hydrological years is required to aid in establishing better watershed-management practices.Several models and indices have been developed for agro-watersheds in temperate climates in Europe, North America and Australia (e.g., [7,[15][16][17][18]), but little work has been done in agro-watersheds in Mediterranean climates which are characterized by a few rainfall-runoff events.

The Concept of P Index
Areas with high potential for P transport not only contain high initial P concentrations but are also characterized by specific land uses and a physiographic nature that easily activates runoff [10].Accurate identification of these areas may enable a change in land use or efficient implementation of remedial options that might decrease P transport from the farmland to waterways.Several models or computational tools have been constructed based on the P index [10,15,[17][18][19], a semi-quantitative tool that includes the various factors influencing the release and transport of P from the terrestrial zone to waterways.The advantage of the P index as a management tool derives from its relative simplicity and the use of existing data sets.The index was first introduced by Lemunyon and Gilbert [20] and was primarily targeted to managing P fertilization in agricultural fields coupled with mapping of their potential P contribution.Later, the index was extended to small agro-watersheds [10] that varied between 0.09 km 2 and 4.5 km 2 [19,21].Heathwaite et al. [15] indicated the need to expand the index to larger watersheds and to test the index on different time scales and during runoff events.
The P index is constructed by weighting the various attributes and factors affecting the amount and transport mechanisms from soils to waterways.In general, each of the attributes is divided into several levels according to their relative contribution to P source and the factors are assigned a weighted coefficient according to their transport potential.Some attributes and factors are specific to site and/or climate; for example, in Nordic areas a freezing attribute was added to account for residual frozen plants in the field [19].Delaune et al. [22] developed an index for pastures in Arkansas, USA, where a multiplication parameter was used to include all of the source, transport and rain attributes.The Natural Resources Conservation Service [23] included soil erosion, runoff, soil P and the amount and method of P addition to the soils.Gburek et al. [10] studied a Pennsylvania watershed and added a component to the index that attested for the hydrological cycle and the distance from the stream to a given area.Heathwaite et al. [15] studied two watersheds in England and devised a scheme that included three hierarchical levels-source layer, transport layer then connectivity layer-which describe the link between the sources and receiving streams.Buczko and Kuchenbuch [17] divided the P index into three computational schemes: the first was an additive approach summing all of the attributes times a weighted coefficient, the second was a multiplication of the source attributes by transport factors, and the third was a mixed approach based on the two previous schemes.The output of all P indices was a map that at the minimum depicted the spatial extents of the CSAs, while the more elaborate computational index yielded a map depicting potential annual P transport from the terrestrial areas of the watershed to the streams.
The main objective of the current work was to develop a P-index model suited to the Eastern Mediterranean climate and the specific land-use practices dictated by the extremely dry summer and wet winter, with most of the precipitation events occurring in the form of rain during a few intense rainfall events [24].In contrast to the above indices, the current work focused on flooding events following our working hypothesis that most P transport in Eastern Mediterranean watersheds occurs during a few extreme events.

Study Area
To test the above hypothesis, we selected a study area within the catchment of Hermon tributary (185 km 2 ) that flows to the Jordan River, which empties into Lake Kinneret.This lake provides about 25% of the state of Israel's potable water and it is known to be P-limited, with periodic cyanobacterial bloom following a decline in N:P ratio [25,26].We assume that the Hermon catchment has high potential for P transport due to high levels of soil P, specific land use (e.g., cattle grazing), low infiltration rates, steep terrain, erodible soils and grazing near stream beds.Flow in the upper reaches of Hermon stream is mostly generated by surface runoff during extreme rain events, whereas the perennial flow in the lower section of Hermon stream is mostly regional base flow.The actual work was limited to two sub-catchments (Figure 1) extending from 1200 m.a.s.l. to 120 m.a.s.l. with average annual precipitation of 600 mm in the lower section and 900 mm at the higher elevation.The number of sampling points along the stream (10 sites) was determined according to the drainage network as well as accessibility to the sampling vehicle.The number of soil samples taken was dictated by soil type, land use and accessibility.The soils of the study area were classified as Chromic Vertisols, Rodexeralfs, Haploxerolls if they exhibited mollicepipedon, and Haploxeralfs [27].

Soil Sampling
We sampled 46 sites with unique georeferencing measured by GPS (E-Trex, Garmin International Inc., Olathe, KS, USA) (Figure 1c).Soil samples were taken from the top layer (0-10 cm) using a composite of four to five sub-samples.The samples were analyzed for bulk density, Olsen-P, water-extractable P (WEP), and texture using standard soil methods [28].The vertical hydraulic conductivity was determined in the field using a mini disk infiltrometer (Decagon Devices Inc., a.

b. c.
Pullman, Washington, DC, USA) [29].To gain a better understanding of the spatial heterogeneity of the soil P, an additional 90 soil samples were collected across the entire Hermon stream watershed and analyzed for Olsen-P and WEP [30].The level of grazing pressure was determined in the field at each of the 46 sites by counting cow droppings along a 50 m × 2 m strip following a modified method discussed by Wood [31] and Madhusudan [32].

Stream Water Analyses
The results of the P-index calculations were tested against stream P concentrations measured at the outlet of each sub-catchment during flooding events (Figure 1b).Water samples were collected in 10 major sub-watershed confluences in the course of four major rainfall events during the 3 years of monitoring that actually generated runoff.These four sampling campaigns yielded 34 water samples with unique spatiotemporal referencing (Table 1).The water discharge was measured with a hand-held Global Water flow probe (model FP-201, Gold River, CA, USA).The relatively small number of events reflected the long period of drought that affected the Jordan River watershed from 2005 to 2011.The water samples were analyzed for total and dissolved P using the ascorbic acid method [30], total suspended solids (TSS) using the gravimetric method [33], and major anions (Cl, NO 3 and SO 4 ) using a Dionex Ion Chromatograph DX-600 equipped with IonPac AS14A anion-exchange column.

Table 1.
Measured discharge values at the 10 sampling stations along the Hermon Stream and its ephemeral tributaries during four hydrological years.

Statistical Analyses
Statistical analyses were carried out using SPSS version 19.The differences in soil P across the study site were tested by ANOVA general model and Kruskal-Wallis test.The coefficient of determination between the P index and P concentrations per event was evaluated by linear regression, while the relationships among the various water constituents were analyzed by Pearson correlation.Spatial statistics within sub-catchments were calculated by Spatial Analyst extension using ArcView Version 9.3.1 (Esri Inc., Redlands, CA, USA).The P index values were aggregated for each subcatchment using spatial mean function.

Arkansas Model
In the current work, we tested several methods of computation and modified existing models to accommodate the conditions of Eastern Mediterranean climate and land-use practices.First, we tested the Arkansas model [22] which was originally designed for pastures and might therefore provide a good solution for our land-use applications (94% being used for cattle grazing, see Figure 1c) with minimal modifications.The values of the model's attributes and a list of the modifications installed in the model to fit the climate and land-use conditions of the study area are summarized in Table 2. Soil-erosion classes were determined according to the RUSLE equation (see below) and the runoff categories were assigned values according to the slope angle and vertical hydraulic conductivity following the soil survey manual [34].For further details of the computation of the P index in the Arkansas model see DeLaune et al. [22].

Revised Universal Soil Loss Equation (RUSLE)
The second model to be tested was the RUSLE, devised by Renard et al. [35].The RUSLE equation was implemented in several P-index models using categorized attributes as demonstrated by Mallarino et al. [36] in Iowa, Good et al. [18] in Wisconsin, and Weld et al. [37] in Pennsylvania.The RUSLE equation was calculated with ArcView 9.3.1 using pixel size of 25 m × 25 m: where A is the amount of erosion per areal unit; R is rainfall erosivity; K is soil erodibility; LS is slope length and steepness; C is crop and management; and P is support method.The R and K values for the Hermon watershed were adapted from Morin et al. [38], LS was derived from digital elevation model (DEM) analysis outlined by Simms et al. [39] and Mitasova and Mitas [40], the C value was taken from the literature [41][42][43] and the P coefficient was assigned a value of 1 since there is no infrastructure in the study area to minimize erosion.For the sake of brevity, the detailed calculation of this model is not presented.

The Runoff Curve Number (RCN) Model
The RCN model was used because its parameters were deemed suitable for the conditions of the study area.The RCN model has been used in Iowa [36] and Canada [44], among other sites.The runoff from a given computational cell within the catchment is computed from a rainfall-runoff equation using ArcCN-Runoff script imported to ArcView 9.3.1.The runoff is computed using the following equation: where Q is runoff; P is precipitation depth; and S is the potential maximum retention of water by soil.Detailed calculation of the ArcCN-Runoff is discussed by Zhan and Huang [45].Within the context of the P loss estimation from soils to waterways the RULSE and the RCN are the transport modules of the P index calculations.

The Hermon Model
The main impetus for the construction of the Hermon model was the partial failure of the above-cited models to accommodate Eastern Mediterranean climate conditions, the lack of some basic data sets, and coarse resolution of existing soil data.The Hermon model was constructed to use actual measurements with little or no categorization.In doing so, we avoided the subjective decision of levels and of setting up boundaries that are not based on theory but on practitioners' decisions.The model's input includes soil P and grazing pressure as source attributes, cumulative rainfall per event (i.e., on average 3 to 4 days) across the catchment, rainfall depth by event, terrain angle, texture, vertical hydraulic conductivity (K) and distance from the streambed as transfer factors (Table 3).The model run begins with an evaluation of the P source, followed by transport computation and concluding with an estimation of potential P transport from a given areal 25 m × 25 m cell.A flow diagram depicted in

Results
The P concentrations according to land-use categories are depicted in Figure 3.The P content in the cattle feeding and resting yards was significantly higher than in the other land-use categories (p < 0.05).The concentrations of dissolved P in the stream water ranged from 0.1 mg•L −1 to 1.6 mg•L −1 while the total P ranged from 0.1 mg•L −1 to 2.2 mg•L −1 .The total P concentrations in the Hermon stream during baseflow conditions are usually lower than 0.1 mg•L −1 , while the soluble P concentrations are less than 0.05 mg•L −1 .The electrical conductance during the runoff events ranged from 160 µS•cm −1 to 712 µS•cm −1 compared to a winter season average of 333 µS•cm −1 at the Hermon spring.In accordance with our working hypothesis, the highest electrolyte concentrations were observed in sampling points adjacent to the cattle feeding and resting areas.For example, the average Cl and SO 4 concentrations in water samples near the feeding areas were 28.5 mg•L −1 and 16 mg•L −1 , respectively, compared with 14 mg•L −1 and 10 mg•L −1 at the other sampling points.
The P-index results from the Arkansas, RUSLE and RCN models are depicted in Figure 4.The darker hues on the maps represent higher P-index values, i.e., a higher likelihood of being a source area.All three models' results suggested that areas adjacent to the streams exhibit a higher potential P contribution.Linear regression analysis showed that the coefficients of determination between the P-index values of the Arkansas and RUSLE models and the measured P concentrations in the streams were rather low (r 2 < 0.32 for both dissolved P and total P), while the coefficients of determination between the RCN model results and the P concentrations in the stream were somewhat higher for dissolved and total P (0.53 and 0.45, respectively).The Hermon model simulations showed that most of the contributing P areas are adjacent to the streams, amounting to less than 15% of the total catchment area (Figure 5).These model results held for both Olsen-P and WEP extracts.The Hermon model results were highly correlated with the measured dissolved-P concentrations in the stream (Table 4).On the other hand, the correlation between the P-index results and the measured total P was less consistent, especially during the rainfall-runoff event of January 2010 where a low r 2 value of 0.14 was computed.

Discussion
The relatively high concentrations of dissolved and total P in the Hermon agro-catchment result from high grazing pressure coupled with numerous cattle feeding and resting spots located near the water source.Similar findings have been reported by McDowell [11] and McDowell and Srinivasan [21] who indicated the urgent need to steer cattle away from streams.Other factors that might create a high potential for P transport from the feeding zones are soil compaction and reduced infiltration, both of which promote runoff [16,46,47].In the Hermon catchment, most of the feeding areas are adjacent to the stream beds, and they therefore receive larger volumes of runoff from contributing areas upslope, resulting in considerable erosional power which mobilizes the P-enriched top layer into the stream.
In the present work we opted to analyze P concentrations in the streams rather than P loads because European and USEPA recommended thresholds of P in the environment are given in concentrations.Moreover, P loads in the streams are strongly affected by autocorrelation which makes statistical analysis difficult.Most of the streams in pasture land of Mediterranean climate are ephemeral in nature which exacerbates the issue of autocorrelation during flow occurrences which are normally happening in response to short burst of rain events in winter time.Hence, in this type of climate the validation analysis of the P index results can only be performed for well-defined flow events rather than annual mean P concentrations or P loads.
The high P concentrations in the stream are about 24 times higher than the total P levels of 0.02 mg•L −1 recommended by the USEPA [48].This threshold is based on the 25th percentiles value for Xeric West Eco-region III while the 25th percentile of the current study was 0.48 mg•L −1 .The high P concentrations in the Hermon tributary during the rainfall-runoff events were also significantly higher than the maximum allowable risk levels in agricultural effluents (0.15 mg•L −1 ) as stated by the European Water Framework Directive-WFD.These results clearly indicate that the cattle-grazing practices in the Hermon watershed are unacceptable in terms of water-quality parameters and more stringent regulations should be enforced to minimize runoff from the feeding and resting hotspots.
The modeling results of the Arkansas, RUSLE and RCN were not correlated with measured P concentrations in the Hermon stream during the rainfall-runoff events.A poor correlation between computed soil-P index and P concentrations in the stream water has also been reported in temperate agro-catchments using these models.For example, Bundy et al. [49] studied an agro-catchment in Wisconsin with a model that was constructed for Pennsylvania by Weld et al. [50] and found a low coefficient of determination (r 2 = 0.14).Once they adjusted the model to the site specific factors affecting the P losses in Wisconsin, a significantly higher correlation was reported (r 2 = 0.79).Bundy et al. [49] explained the differences in model results as the outcome of using different erosional coefficients in the RUSLE equation which is embedded in these models.On the other hand, Harmel et al. [51] studied an agro-catchment in Texas and showed that the main source of error in the P-index calculations resulted from their inability to identify the various erosional levels.This type of error was minimized in the model developed by Harmel et al. [51] because they opted to use measured data rather than estimating the erosional levels.
The Hermon model used in the current work employed data measured in situ such as texture, and hydraulic conductivity rather than values imported from the literature.Moreover, we avoided the subjective categorization of the continuous variables of the model as practiced in previous studies.The high correlation between the P values computed by the Hermon model and the measured P concentrations in the stream strongly support this modeling approach.Harmel et al. [51] compared the Arkansas, Iowa and Texas models and found a high correlation between the index results and P concentrations and loads using the multiplicative Iowa model (r 2 = 0.91), a moderate correlation with the additive Texas model (r 2 = 0.5) and a low correlation with the Arkansas model (r 2 = 0.3).Their results were in agreement with the current work which emphasizes the importance of using measured site-specific factors affecting erosion rather than general values from the literature.
The correlation between the computed P-index values and P concentrations in the stream during the rainfall-runoff event of January 2010 was rather poor compared with the other events monitored during the 3-year study.The reason for the low predictive power of the model in this event was the high resuspension rate of stream sediments.This particular rainfall event generated the first significant runoff of the hydrological year which came after a long hiatus in runoff events.Huisman and Karthikeyan [52] found that small runoff events tend to deposit suspended-P material on the stream bed.This material is then remobilized into the water column during the turbulent flow characteristic of major runoff events (see Table 1), consequently increasing the total P concentration in the stream.The loss of predictive power during the first major rainfall-runoff event of the hydrological year attests to the importance of correctly structuring the model.Vadas et al. [53] formulated a procedure for optimal model structuring to problems similar to P-index computation.The procedure begins with a perceptual model that describes the general qualitative understanding of the problem; this qualitative formulation is then translated to a conceptual model which consists of a set of equations and continues to a numerical solution using a well-designed algorithm.A flawed, inaccurate or incomplete perceptual model will result in an unreliable prediction.In recent years, a framework for diffuse pollution known as the transfer continuum [54] has been widely used as a perceptual model for P mobilization and transport.In fact, most of the P-index models utilize the transfer continuum to produce a set of equations that aim to identify CSAs and can be viewed as a conceptual model ( [10,15,17,19,49,51] among many others).The current study clearly demonstrates that adapting the perceptual transfer continuum model to a Mediterranean catchment fails to consider P resuspension in streams during high-flow events.Hence, the modeling procedure has to be revisited and a new set of equations added to account for the P-transport mechanism in the streams, which is highly dependent on rainfall duration, intensity, frequency and soil antecedent moisture that can generate considerable runoff at the beginning of the hydrological year.
The main sources of uncertainty in the Hermon model are the omission of temporal variability of factors affecting P losses within and between the hydrological years.These sources of uncertainty include groundcover, antecedent soil moisture and unaccounted changes in land use generated by farming activities.Future work should consider the inclusion of remote sensing information on groundcover and soil antecedent moisture variability within the spatially identified critical source areas.These measurements undoubtedly would improve the calibration and verification of the P model.

Conclusions
Several leading P-index models were tested toward our main objective of developing a P-index model suitable for the Eastern Mediterranean climate.None of the published models performed well under the climatic conditions of extremely dry summers and wet winters, where most of the precipitation events occur in the form of rain during a few intense rainfall events.We developed a new P index termed Hermon model which was adjusted to the rainfall-runoff event occurrences.In general, extremely high correlations were recorded between the P-index values and the P concentrations in the stream.One notable exception was recorded in the correlation matrix which indicated that the perceptual model that underlies this computational approach requires some modifications to address the resuspension of P from the stream bed during the typical early runoff events in Eastern Mediterranean catchments.The results also suggested that the current cattle-grazing practices in the study area are unacceptable in terms of water-quality parameters and more stringent regulations should be enforced to minimize the runoff from the feeding and resting hotspots.

Figure 1 .
Figure 1.(a) The general location of Hermon stream; (b) the locations of the studied sub-catchments; and (c) the current land use and stream sampling sites.

Figure 2
Figure 2 describes the sequence of the calculations along with the algorithm summarized in the equations part of the figure.

Figure 3 .
Figure 3. Distribution of (a) water-extractable P (WEP); and (b) Olsen-P in the four selected land-use areas of Hermon catchment.

Table 4 .
Regression analysis results between P concentrations in the streams and P-index values computed by the Hermon model.

Figure 5 .
Figure 5.The Hermon model P-index map.The model emphasizes the importance of the areas located close to the stream beds, while the upper part of the basin is less likely to contribute P to the waterways.

Table 2 .
Values of the Arkansas model's attributes.The values were modified to accommodate the Eastern Mediterranean conditions as follows: (1) soil test was changed to water-extractable P; (2) grazing density was used as an indicator for added P; (3) changes in precipitation depth; and (4) no best management option was included.
I WSP-water soluble P.

Table 3 .
Attributes used to compute the P index with Hermon Model.The soil-P represents average values in mg/kg.
Figure 2. Flow diagram of P-index calculations using ArcView model builder.