MODELING OF PHOSPHORUS LOSSES BY WATER EROSION

Phosphorus losses in runoff water and eroded sediment may increase the risk of eutrophication. Erosion prediction models can be used to assess and quantify nutrient losses and transport in different soil management scenarios. This study aimed to assess the performance of models to estimate the losses of available phosphorus in eroded sediments and assess their spatial variability pattern. The experiment was installed on a eutrophic Red Ultisol located in Tabapuã, northwestern São Paulo State, São Domingos Stream Watershed. Aiming at estimating and validating the models, measurements were made from 2007 to 2015 at 17 observation points (slopes), determining soil phosphorus content (Psoil), eroded sediments (Psed), and enrichment rate (ERsed). The models applied to estimate the enrichment rates of phosphorus in the eroded sediment and runoff are efficient when in use with the predictions of the model Water Erosion Prediction Project (WEPP). The enrichment rates obtained presented a Nash-Sutcliffe (NS) coefficient close to 1. Losses of soluble phosphorus in runoff above 0.02 mg L, a critical value for eutrophication, can be obtained in 81% of the watershed area with a probability higher than 75%.


INTRODUCTION
Water erosion is one of the main problems of soils and lands worldwide and in Brazil since it involves detachment, transport, and deposition of soil particles, also transporting adsorbed organic and chemical elements to aquatic environments.
In this sense, phosphorus stands out among the elements adsorbed in the sediments.Excess nutrients may cause environmental imbalance since both runoff and sediments act as natural and important vectors for the transfer of these elements to the aquatic biota, in addition to causing sedimentation and eutrophication, which compromise the watercourses (Chakraborty et al., 2014;Bhadha et al., 2016).
Considering the diversity of research methods to assess the occurrences and consequences of the water erosion process, prediction models are essential since they can be used for soil loss and sediment transport studies in different soil management scenarios due to the union of knowledge and data generated by the physical models.With these models, conservationists will be led to the choice of conservation practices able to fit the needs of a specific area (Pieri et al., 2007;Flanagan et al., 2013).
An erosion prediction model is represented by a set of equations consisting of factors such as climate, soil, topography, soil management, soil loss, and sediment (Flanagan et al., 2012).Models such as the Water Erosion Prediction Project (WEPP) have provided good results in predicting erosion rates (Amorim et al, 2010;Mahmoodabadi & Cerdá, 2013).WEPP is widely used by researchers in temperate regions (Foltz et al., 2011;Pieri et al., 2014).For this reason, verifying the applicability of this model to the Brazilian edaphoclimatic conditions is essential.This verification has the purpose of assisting in agricultural planning related to soil and water conservation, respecting the limits established by legislation.
In Brazil, in most studies of water erosion, WEPP model has traditionally been used to monitor and represent erosion spatially and temporally, with details of the processes of transport and deposition of sediments in slopes or small watersheds.In this context, the scientific and technological challenge of this study is related to the adaptation and application of WEPP in modeling the transport of phosphorus.This will allow verifying possible impacts from land use and management changes in agricultural soils on the process of eutrophication of water bodies by the presence of phosphorus.
Engenharia Agrícola, Jaboticabal, v.38, n.1, p.149-157, jan./feb. 2018 Accordingly, the aim of this study was to assess the performance of the WEPP model to estimate the losses of available phosphorus in eroded sediments and assess their spatial variability pattern.

MATERIAL AND METHODS
The experimental area was located in Tabapuã, northwestern São Paulo State (21°05′57″ S and 49°01′02″ W) (Figure 1a).Regional climate is a C2dA′a′ type according to Thornthwaite (1948), which stands for a megathermal climate, with an average annual precipitation of 1,318 mm.
Regional primary vegetation was classified as a seasonal rain forest and Cerrado.The area has 220 ha and presents a history of more than 20 consecutive years of sugarcane (Saccharum spp.) cultivation, with more than 10 years under the mechanized harvesting system without burning.Sanchez et al. (2009) classified the soil of the area as a eutrophic Red Ultisol, according to Embrapa (2013) criteria.Two landforms, a concave and other convex, were identified in the area (Figure 1b).Particle size distribution and mineralogical composition are shown in Table 1, according to Barbieri et al. (2009) and Camargo et al. (2013a;2013b).Erosion and sediment production were estimated by WEPP on 89 slopes selected by Ibiapina (2015) from experimental results obtained by Martins Filho (2007).The 89 slopes included a set of 626 points spaced every 50 m.The values of phosphorus concentration in the soil (Psoil) and in the eroded sediments (Psed) were obtained from 17 observation points (slopes).The enrichment rate (ER) was obtained from the relation between Psed and Psoil.Five types of input information were used in the WEPP model: (i) climatic, (ii) hydrological, (iii) physiology of sugarcane crops, (iv) soil type, and (v) erosion/deposition data.Climatic information was obtained stochastically by using the parallel model CLIGEN to generate data from observations of precipitation, temperature, dew point, and solar radiation of the meteorological station of Catanduva, São Paulo State, where data are complete and similar to Tabapuã.
Hydrological information was generated from the characteristics of evapotranspiration, topography, and infiltration (Figures 1b, 1c, and 1d).Soil information (texture, density, organic matter content, and erodibility) were obtained from Martins Filho (2007) and Sanchez et al. (2009).Particle size classification on input information on soil followed the USDA (2014) classification criteria.
Physiological characteristics of sugarcane crop were those used by Ibiapina (2015) in a WEEP model.In the WEPP scenarios, soil cover was equivalent to 90%.For the input information on soil management, we used the values described by Pieri et al. (2007).
The composition of the eroded sediment was classified as follows: 1) primary fractions (<0.002 mm), sand (0.06 to 2 mm), and silt (0.002 to 0.06 mm); 2) small aggregates consisting of clay, silt, and organic matter; large aggregates constituted by the three primary fractions and organic matter (Elliot et al., 2015).
For the fractionation of the eroded material, the approach proposed by Foster et al. (1985) and incorporated in the WEPP was used.This fractionation allowed obtaining the enrichment rate WER (Equation 4) for validation tests in 17 slopes from those 89 established in the study area (Figure 1b).
The estimate of phosphorus transported with eroded sediments was obtained by adopting two algorithms.In the first algorithm, eqs (1) and (2) were used: where, Psed is the amount of phosphorus in the sediment (kg ha −1 ); Psoil is the content of phosphorus in the soil (kg kg −1 ); ER is the phosphorus enrichment rate, and Sed is the amount of eroded sediment (kg ha −1 ).
The values initially used for a and b were 2.0 and −0.2, respectively (Perez-Bidegain et al., 2010).Both a and b values were calibrated with the ER values obtained on seven slopes and, subsequently, the performance of the Engenharia Agrícola, Jaboticabal, v.38, n.1, p.149-157, jan./feb. 2018 resulting equation was tested by using ER from 10 slopes.
In the second algorithm, the amount of phosphorus transported with eroded sediments was estimated according to eqs (3) and ( 4): WER is the phosphorus enrichment rate computed by WEPP; SSAsed is the specific surface area of the sediment (m 2 g −1 ), and SSAsoil is the specific surface area of the original soil (m 2 g −1 ).
The performance of the studied models was assessed by using different criteria, which included the Nash-Sutcliffe (NS) coefficient, root mean square error (RMSE), values of mean of errors or mean differences (MD), and concordance index (d), according to eqs ( 5), ( 6), ( 7) and ( 8 The maximum value for the NS coefficient is 1.MD and NS can be negative.If the NS coefficient is lower than zero, the model predicts values worse than the observed average.When MD presented a (+) or (−) sign, we assumed that the predicted values, on average, overestimated or underestimated the observed values, respectively.The index d is delimited between 0 and 1.
Geostatistical analysis was performed following the methodology described by Vieira et al. (1983) and Robertson (1998).In order to fit the semivariogram, 366 points were randomly chosen in the space, as performed by Gertner (2003).The other 260 points were used to compare validation parameters between the ordinary kriging (KO) and sequential Gaussian simulation (SGS) methods.

RESULTS AND DISCUSSION
The average amount of eroded sediment estimated with the WEPP model was 3097.63 kg ha −1 , considering a mechanized harvesting without burning.This amount was lower than the annual average estimated by Sanchez et al. ( 2009) by using a USLE for an Ultisol in Tabapuã, São Paulo State, reaching a value of 8260.00 kg ha −1 in a sugarcane field harvested manually after burning.Martins Filho et al. (2009) obtained average losses of 2836.67 kg ha −1 in experimental plots under simulated rainfall in an Ultisol area under the mechanized harvesting system without burning.This indicates a possible overestimation by WEPP of this last observed value.
The average contents of phosphorus available in both soil and eroded sediment (Psoil and Psed) correspond to the class medium (16-40 mg dm −3 ) for sugarcane cultivation (Van Raij et al., 1998).These contents are similar to those obtained in studies carried out on soils under sugarcane crop cultivation (Camargo et al., 2013a;Oliveira et al., 2013).On average, eroded sediments were 1.51 times more phosphorus-enriched than the original soil.Silva et al. (2012) observed average losses of available phosphorus in an Oxisol under sugarcane cultivation of 26 mg dm −3 for an enrichment rate of 1.28.These results are lower when compared to those obtained in our study, but consistent with the more erodible nature of the analyzed Ultisol.Martins Filho et al. (2009) observed enrichment rates of phosphorus-eroded sediment from 2.7 to 3.8 in experimental plots of 35 m 2 , being higher than those assessed in our study (Table 2).
Engenharia Agrícola,Jaboticabal,v.38 Table 2 shows that the values of the enrichment rate WER, obtained by means of WEPP and the adjusted function ln(ER) = 2.682 − 0.278 ln(Sed) (R 2 = 0.69, p<0.02), were, on average, significantly higher than those obtained with the equation ln(ER) = 2 − 0.2 ln(Sed), as presented by Perez-Bidegain et al. (2010).In order to better assess and understand the obtained results, as well as the reliability of the model estimates, we used statistical indices such as the mean difference (MD, Table 3), which shows the trend of error or vice of the models.Nash-Sutcliffe coefficient values indicate the fit of the simulated data to those observed in relation to the 1:1 straight line, which may vary from −∞ to 1. NS is associated with the estimated efficiency of phosphorus enrichment in the eroded sediment.As described by Lopes et al. (2016), NS ≥ 0.75 represents a model considered as adequate or good, 0.36 < NS < 0.75 represents a model considered as acceptable, and NS < 0 represents a model in which the data series is best represented by the mean use.In this sense, the three tested models were classified as adequate or good.
In the analysis of the root mean square error (RMSE), model 1 presented a lower value when compared to models 2 and 3 (Table 3).RMSE values of 0 are considered equal.
Values of MD > 0 indicate an overestimation of phosphorus values by the model when compared to the observed or measured values, as in models 1 and 2. Model 3 tends to underestimate the observed values.Because MD is a measure of accuracy, which is determined by Equation ( 7), if its signal is negative, the estimated value will be lower than the observed, which will indicate a tendency of the fitted model to underestimate true values.
The concordance index (d) of Willmott (1982), which expresses accuracy between estimated and observed values, is satisfactory (d > 0.75) for the three models, according to the criteria adopted by Lopes et al. (2016).The models provided a high degree of accuracy since their concordance indices (d) were close to unity.
Our results are in accordance with those conducted by Perez-Bidegain et al. ( 2010) and Elliot et al. (2015), showing the possibility of incorporating into WEPP an algorithm, or routine to predict phosphorus derived from slopes with eroded sediment.Thus, the spatial distribution of available phosphorus losses in the eroded sediment was estimated for the watershed area composed by 89 slopes.For this, the models presented in Table 3 to estimate the phosphorus enrichment rate of the eroded sediment were used.
The enrichment rates WER did not differ significantly between convex and concave slopes (Table 4).However, enrichment rates on concave slopes were significantly higher when compared to those convex, when obtained with the methods ( 2) and (3).
Engenharia Agrícola, Jaboticabal, v.38, n.1, p.149-157, jan./feb. 2018 When the phosphorus originally available in the soil was analyzed, no significant difference was found between the convex and concave landforms.In the sediment and runoff water, the average concentrations of available phosphorus were significantly higher in concave slopes than in the convex slopes (Table 4).This result might be attributed to the higher production and sediment deposition on the concave slopes.In the processes of water erosion, P losses occur due to the adsorption of this element on the sediment surfaces or in suspension in the water flow.
All averages of soil erosion variables such as net erosion rate (A), sediment concentration in the runoff (SC), runoff volume (RV), and production and deposition of sediments were significantly higher on the concave slopes when compared to those convex (Table 4).These results contradict those of Sanchez et al. (2009), which were obtained in the same area of our study.Sanchez et al. (2009) estimated, by using the Universal Soil Loss Equation (USLE), higher soil losses due to erosion for convex landforms in relation to those concave.In fact, USLE does not consider the sediment deposition process, in addition to not estimating the furrow erosion in a more advanced phase, which may explain the results found by Sanchez et al. (2009).Another limitation in using USLE is the calculation of the relief factor (LS) in convex ramps for use in watersheds.In this case, the larger the ramp length is, the greater the soil loss.Moreover, Ibiapina (2015) demonstrated that the WEPP model provided a greater precision and accuracy in estimating the erosion values in relation to USLE.WER and ERenrichment rates; P(1), P(2), and P(3)concentrations of available P in the sediment and runoff for enrichment rates (WER and ER) obtained with: (1) WER of the WEPP model, (2) ln(ER) = 2.682 − 0.278 ln(Sed), and (3) ln(ER) = 2 − 0.2 ln(Sed), respectively; Aerosion rate; SCsediment concentration in the runoff; RV − runoff volume.Means followed by the same letter in the row do not differ from each other by the t-test (p<0.05).Camargo et al. (2013b) assessed the mineralogy of the Ultisol used in this study, in relief landforms, and found the highest contents of hematite and goethite in the convex area, with a predominance of hematite in both landforms (convex and concave).In the convex area, hematite presented the smallest average diameter of crystal, with a predominance of kaolinite in the concave area.
The presence of kaolinite implies the development of a block-type macrostructure, as observed by Ferreira et al. (1999), originating a soil with higher density, a higher proportion of small pores, and lower permeability.Therefore, erodibility should be higher in kaolinitic soils.This could explain the highest values of A (Figure 2), SC, RV, and a production and deposition of sediments significantly higher in the concave slopes when compared to those convex.For the Ultisol assessed in this study, Barbosa (2015) demonstrated that furrow erodibility values (Kr) in the convex landform are lower when compared to the Engenharia Agrícola, Jaboticabal, v.38, n.1, p.149-157, jan./feb. 2018 concave landform, with values of 0.0053 kg N −1 s −1 and 0.0146 kg N −1 s −1 , respectively.The critical shear stress values are 2.1887 N m −2 in the convex landform and 2.7260 N m −2 in the concave landform.These results are consistent with the erosion patterns established by the WEPP model in the studied watershed and obtained by Martins Filho (2007), in addition to corroborating the results of Camargo et al. (2013b).Local mineralogy has an influence on soil attributes and local hydrology, which affects the degree of local erosion.Thus, this will also influence phosphorus adsorption to the eroded sediment.
As showed in the results, the highest phosphorus losses occurred in concave slopes.According to Camargo et al. (2013a;2013b), the convex area showed the highest values of weighted mean diameter, geometric mean diameter, aggregates >2 mm, aggregates of 2 to 1 mm, total pore volume, and soil water content; as well as the lowest values of soil density and resistance to penetration at the two studied depths (0.0-0.2 and 0.2-0.4m) when compare to the concave area.Camargo et al. (2013b) also observed an expressive number of soil attributes with a greater spatial variability in the concave area.
A relevant aspect is the influence of clay fraction mineralogy on P adsorption, as pointed out by Camargo et al. (2013a).Phenomena such as the adsorption of elements in the soil are strongly related to mineral crystallinity.Thus, phosphorus tends to be adsorbed by goethite (Gt), which has a low crystallinity and high SSA, as well as by gibbsite (Gb).This allows the use of the technique described in this study as an alternative to management practices in order to minimize the effects of phosphorus adsorption in the soil, for example.
Among the three models tested to estimate the enrichment rate (ER or WER), no difference was observed in the spatial pattern of losses of the available P in the watershed area (Figure 3).Enrichment rate was estimated by WEPP, i.e.WER (Figure 3a), as well as for empirical models (Figures 3b and 3c).
The annual average concentrations of available or soluble phosphorus in the runoff, derived from 17 slopes in the watershed, ranged from 0.059 to 0.067 mg L −1 (Table 5).These values are worrisome because they are three times higher than the critical limit for the total phosphorus present in water (0.020 mg L −1 ).Therefore, the use of models that estimate not only the surface runoff and the sediment derived from the slopes, but also assess P losses in the sediment and runoff water is essential, as observed by Elliot et al. (2015).Engenharia Agrícola,Jaboticabal,v.38,n.1,The spatial distribution pattern of concentrations of P > 0.020 mg L −1 was similar between models 1, 2, and 3 (Figures 4a, 4b, and 4c).For these models, in the most critical regions (p>0.75), the probability of P > 0.020 mg L −1 occurs in more than 81% of the watershed area (Table 6).For probabilities (p) <0.25, 0.25 to 0.50, and 0.50 to 0.75, the percentage of representative areas occupied by these classes in the watershed are, on average, 1.7, 6.6, and 10%, respectively.
The eutrophication of water bodies by phosphorus is considered a great concern, as shown by Wang et al. (2011).According to the Resolution of CONAMA No. 357/2005, the critical level of total phosphorus in water is 0.020 mg L −1 for lentic environments (Brasil, 2005).Thus, the transfer of phosphorus from the terrestrial system to the aquatic environment is a key factor and occurs mainly in two ways: surface runoff and percolation in the soil profile.Phosphorus forms transferred to aquatic environments can be soluble and particulate.Particulate phosphorus is bound to mineral and organic colloids, characterizing phosphorus as inorganic and organic (Klein & Agne, 2012).Engenharia Agrícola, Jaboticabal, v.38, n.1, p.149-157, jan./feb. 2018 TABLE 6. Percentage of the watershed area by probability interval (p) for the available phosphorus concentration be higher than 0.020 mg L −1 in the runoff.

CONCLUSIONS
The models used to estimate the enrichment rates of the eroded sediment and runoff are efficient when in use with the predictions of the model Water Erosion Prediction Project (WEPP) and present a Nash-Sutcliffe (NS) coefficient close to 1.
The use of spatial distribution of loss of available phosphorus with eroded sediment in the study area allows identifying specific management areas and sites with a greater loss of phosphorus.
The highest losses of phosphorus occur on slopes with concave landforms, reaching values higher than 0.020 mg L −1 , a critical value for eutrophication, and can be obtained in 81% of the watershed area with a probability greater than 75%.

FIGURE 1 .
FIGURE 1. Location of the area close to the Turvo river (a), slope map and 17 observation points for simulated rainfall information collection (b), surface water flow simulation (c), and planialtimetric map with altitude in meters (d).Adapted from Ibiapina (2015).
value estimated by the model, Oi is the observed or measured value, n is the number of observed values, and O is the average of the measured values.

FIGURE 2 .
FIGURE 2. Spatial pattern obtained by ordinary kriging of the estimated net erosion (A) with the WEPP model.

FIGURE 4 .
FIGURE 4. Spatial distribution obtained by sequential Gaussian simulation of the probability of the available phosphorus concentration be greater than 0.020 mg L −1 in the runoff, with values of enrichment rates (WER and ER) obtained as: (a) WER of the WEPP model, (b) ln(ER) = 2.682 − 0.278 ln(Sed), and (c) ln(ER) = 2 − 0.2 ln(Sed).

TABLE 1 .
Particle size distribution and mineralogical composition of soil in the study area.

TABLE 2 .
Amount of eroded sediment, phosphorus (P) concentrations in the soil and sediment, and P enrichment rates for 17 slopes in the study area.Sedamount of eroded sediment estimated by WEPP derived from each slope (kg ha −1 ); Psoilsoil phosphorus content (kg kg −1 ); Psedamount of P in the sediment (kg ha −1 ); ERphosphorus enrichment rate; Nenon-estimated value.(2) ER obtained with the adjusted function ln(ER) = 2.682 − 0.278 ln(Sed) (R 2 = 0.69, p<0.02), obtained with data from slopes highlighted with *. (3) ER obtained with the function ln(ER) = 2 − 0.2 ln(Sed).Means followed by the same letter in the row do not differ from each other by the t-test (p<0.05).

TABLE 3 .
Performance of three methods for estimating the phosphorus enrichment rate in the eroded sediment.

TABLE 4 .
Phosphorus in the soil, enrichment rates, phosphorus in sediment and runoff water, and soil erosion variables as a function of relief landform.

TABLE 5 .
jan./feb.2018 Concentrations of eroded sediment and phosphorus in the soil and runoff water on 17 slopes of the watershed.