Particulate PhozzyLogic Index for policy makers—an index for a more accurate and transparent identification of critical source areas

This study presents an algorithm for the allocation of particulate phosphorus (PP) loads entering surface waters to their sources of origin, which is a basic requirement for the identification of critical PP source areas and in turn a cost-effective implementation of mitigation measures. Furthermore, it conducts a sensitivity analysis determining the impacts of storm drains, discharge frequencies and flow directions on the designation of critical source areas with the help of present-day scenarios for a case study catchment with an area of several hundred square kilometres. The results show that the scenarios at least partially disagree on the identified critical source areas and suggest that especially open furrows at field borders have the potential to cause deviating critical source areas. Finally yet importantly, the Particulate PhozzyLogic Index based on model results as well as fuzzy logic is introduced. Its main feature is to transform the results of diverse scenarios or even models into a final map showing a catchment-wide ranking of the possibility of high PP emissions reaching surface waters for all agricultural fields.


Introduction
An innovative aspect of the European Union Water Framework Directive (European Comission, 2000) is its inclusion of economic principles in the achievement of its objectives (Martin-Ortega and Balana, 2012).One of these principles is cost-effectiveness.Mitigation measures capable of contributing to the achievement of the "good status" have to be analysed and compared in the light of it.The effectiveness and therefore cost-effectiveness of mitigation measures, however, is not uniform and strongly depends on the location of their implementation.
This non-uniformity exists on different scales and ranges from the national to the local scale and beyond.
In the field of diffuse pollution, one can find the concept of critical source areas on one of the largest scales, i. e. often the field scale.Critical source areas are commonly defined as those source areas, which dominantly contribute to the pollution load of surface waters via one of many possible transport pathways (e. g.Heathwaite et al., 2005b;Pionke et al., 2000;Strauss et al., 2007).This concept is especially applicable to sediment-bound pollutants like particulate phosphorus (PP) where usually only a small portion of a watershed is responsible for the majority of phosphorus inputs into its surface waters (e. g.Kovacs et al., 2012;Strauss et al., 2007;White et al., 2009).
Since the effectiveness of mitigation measures is highest on critical source areas, their identification takes a key role in a cost-effective approach to evaluate and implement best management practices.In this process, models and modelling studies are an important factor (Cherry et al., 2008;Sharpley et al., 2015) and tiered approaches have proven useful in limiting efforts required for modelling (Doody et al., 2012;Heathwaite et al., 2005a;Zessner et al., 2014).Heathwaite et al. (2005a), for example, applied a very basic risk screening approach depending only on two variables to all lakes in Great Britain and recommend the application of a semi-distributed indicator tool to previously revealed high-risk catchments.Sharpley et al. (2009) state as a rule of thumb that about 80 % of the phosphorus inputs into the surface waters of a catchment originate from only about 20 % of its area.This relationship is sometimes also known as the 80:20 rule or the Pareto principle.An early approach to identify such hotspot areas is the Phosphorus Index (Lemunyon and Gilbert, 1993), which was originally developed in the USA but has found its way to Europe since then (Heathwaite et al., 2003).By now, there exist several modifications and/or localisations of the Phosphorus Index (Buczko and Kuchenbuch, 2007).A general advantage of all of them is that they are simple, good to communicate and are capable of identifying high-risk regions for a cost-effective implementation of mitigation measures (Cherry et al., 2008).
An important modification comprises the separation of source and transport factors and the incorporation of hydrologic return periods (Gburek et al., 2000).
Although the Phosphorus Index is basically a semi-quantitative tool, especially the Iowa Phosphorus Index tries to approximate biologically available phosphorus loads entering surface waters (Mallarino et al., 2002).Beyond that, with its stronger physical basis, the Swedish Phosphorus Index requires more input data and is not as easy to calculate as other Phosphorus Indices (Buczko and Kuchenbuch, 2007).This leads us to another approach to identify critical source areas: semiempirical/conceptual, spatially distributed soil erosion and phosphorus transport modelling.While this approach usually requires even more input data and higher calculation efforts, it can be expected to provide results, which are more accurate, due to its improved quantitative basis compared to the Phosphorus Index approach.Nonetheless, such models generally represent a good compromise between solely empirical and process based models (Cherry et al., 2008).Respec-tive studies have been carried out utilising, for example, the WaTEM/SEDEM (Onnen et al., 2019;Van Oost et al., 2000;Van Rompaey et al., 2001;Verstraeten et al., 2002) and PhosFate (Kovacs et al., 2008(Kovacs et al., , 2012) ) models in the past.
One major drawback of these modelling studies is, however, that it is not outright clear how to translate model results into catchment-wide risk assessments on field-level.We therefore propose a novel semi-empirical "index" for PP based on the results of spatially distributed models, i. e. PhosFate in our case, as well as fuzzy logic (Particulate PhozzyLogic Index -PPLI), thus trying to overcome this shortage and aiming at reducing vagueness in identifying critical source areas at the same time.Fuzzy logic was developed in the mid-sixties of the 20 th century and is based on fuzzy sets (Zadeh, 1965).Later it was developed further into a theory of possibility (Zadeh, 1978) and Singpurwalla and Booker (2004) yet show a way of incorporating fuzzy sets into probability theory.
Although tempting, we do not further explore the latter concept, as possibility theory and Phosphorus Index surprisingly share some very basic ideas in the first place.We hypothesise that the potential of phosphorus loss of or movement into surface waters from a certain site as it is used by the Phosphorus Index is conceptually related to possibility theory yet lacking its theoretical foundation in mathematics.Hence why it is called an index.Robinson (2003) provides an overview of fuzzy sets and on how they are used in modern geographic information systems (GIS).Incorporating all this into probability theory in turn rather is capable of confusing policy makers facing the task to meet water quality targets in a certain catchment than it does add to dealing with uncertainty when it comes to this.
While fuzzy logic previously has been successfully applied to, for example, river quality (e. g.Lermontov et al., 2009;Liou et al., 2003), groundwater quality (e. g.Rebolledo et al., 2016;Vadiati et al., 2016) and landslide susceptibility (e. g.Champati ray et al., 2007;Pourghasemi et al., 2012;Tien Bui et al., 2012) analyses, to our knowledge, it has not been used in the context of the identification of critical source areas so far.
This study thus tries to assess the possibilities and limitations of combining an enhanced PhosFate model with fuzzy logic for the identification of critical PP source areas.It does so by conducting a sensitivity analysis determining the impacts of selected uncertain input data as well as model parameters on the designation of critical source areas in the form of present-day scenarios for a case study catchment with an area of several hundred square kilometres.
As an important prerequisite for this serves the presented model's enhanced algorithm for allocating the PP emissions actually reaching surface waters to their respective source areas.

Case study catchment
The present study was conducted on the catchment of the river Pram in the north-western part of the Austrian federal state of Upper Austria.According to  Kovacs et al. (2008Kovacs et al. ( , 2012)), the semi-empirical, spatially distributed phosphorus emission and transport model PhosFate was recently extended by Hepp and Zessner (2019) with a module capable of taking storm drains at road embankments into account.It models erosion with the help of the (R)USLE (Renard et al., 1997;Schwertmann et al., 1987;Wischmeier and Smith, 1978) making use of a raster GIS data based single (D8) flow algorithm version of the slope length factor (Desmet and Govers, 1996).The PP emission part of the model then combines the erosion with the PP content of the top soil and a local enrichment ratio of each raster cell in order to calculate local gross PP emissions (Kovacs, 2013).
PP retention in turn is computed via an exponential function of the cell residence time and a mass balance equation including terms for the inflowing PP load, the local gross PP emission and the local as well as the transfer PP retention.The hydraulic radius, among others, is a requirement for the calculation of the cell residence time and itself involves model parameters related to discharge frequency (Kovacs, 2013).With the channel as well as overland PP deposition rates and the PP transfer coefficient of the storm drains model extension, PhosFate features three potential calibration parameters.

Revised allocation algorithm
One major output of PhosFate is the calculation of PP cell loads via an allocation algorithm.The PP cell loads describe the amount of local PP emissions actually reaching surface waters, i. e. the PP cell load of a single cell represents its local PP emission minus the cumulated amounts of retention taking place in all of its downstream cells (Kovacs, 2013).Since the original allocation algorithm of the PhosFate model based on transmission coefficients was developed for larger cell sizes (e. g. 100 × 100 m resolution) than used in this study (10 × 10 m resolution), it was necessary to revise this algorithm.
While the original algorithm only required a single top-down computation starting with the lowest and finishing with the highest flow accumulations, the revised algorithm builds on the existing top-down computation resulting in the PP retentions as well as transports and adds an additional bottom-up computation starting with the highest and finishing with the lowest flow accumulations to it.The latter results in the PP cell loads for which the local net PP emissions calculated as local PP emissions minus PP retentions constitute the upper PP cell load limits.This computation in turn involves the calculation of PP cell transfers.These are the cumulated amounts of PP entering a cell from upstream cells, which are transferred through the cell and actually reach a surface water (again minus the cumulated amounts of retention taking place in all of its downstream cells).The complete bottom-up computation thus consists of cycles of the following steps, which, however, are executed for overland cells only: 7. Weighted apportioning of the cell load carry-over to the inflowing cells utilising their PP transports as weights.

Coupling of the PhosFate with the STREAM model
The STREAM model was designed to simulate overland flow and erosion in agricultural catchments (Cerdan et al., 2002a,b;Le Bissonnais et al., 2005).Yet it does lack the ability to simulate the emission and transport of chemical substances like phosphorus.What makes it appealing from the point of view of phosphorus emission and transport modelling, though, is its capability to simulate a flow network, which accounts for tillage directions among others (Couturier et al., 2013).The effect of tillage direction on flow direction has been thoroughly studied by Souchere et al. (1998); Takken et al. (2001a,b,c) and depends mainly on the angle between tillage direction and topographic slope aspect, topographic slope gradient and surface roughness.
In order to let PhosFate benefit from the last mentioned feature, we coupled both models, i. e. we used version 3.7.1 of the STREAM model to pre-process the flow direction data handed over to PhosFate.For this, we applied a stepwise approach with (i) topographic flow directions only, (ii) the sole enforcement of tillage directions on topographic flow directions for arable land and (iii) the combined enforcement of tillage directions and open furrows at all field borders (sides as well as headlands, which demonstrates an extreme case) on topographic flow directions for arable land.Fortunately, field border data was available from a (geo-)database related to the Integrated Administration and Control System (IACS) of the Common Agricultural Policy (CAP) of the European Union (EU) (Hofer et al., 2014).
As tillage directions are an input parameter to the STREAM model, they first had to be derived.This was accomplished by means of minimum bounding rectangles with smallest areas enclosing each field and the assumption that tillage directions are parallel to the longer two sides of those rectangles.A sample of 176 fields was then used to compare such derived tillage directions to visually determined tillage directions from orthophotos.In only about 16 % of the cases (29 fields) they deviated by more than ±22.5 • .Two more required input parameters are the average surface roughness in tillage direction and the average surface roughness perpendicular to tillage direction.These were globally set to 1 to 2 cm and 2 to 5 cm respectively.Slopes were generally calculated in the direction of the topographic or enforced D8 flow directions.

Present-day case study scenarios
Drained roads and storm drains can be a relevant emission pathway into surface waters in Switzerland (Alder et al., 2015;Bug, 2011;Doppler et al., 2012;Prasuhn, 2011;Remund et al., 2021).Hepp and Zessner (2019) come to the same conclusion for Austria and estimated that in the same case study catchment as examined in this study about 77 % of all fields upstream to roads are artificially drained by one or more storm drains with a 90 % highest posterior density interval (credible interval, sometimes also called the Bayesian confidence interval) ranging from approximately 54 % to almost 100 %.They further performed a plausibility check of the transfer coefficient of the storm drains model extension accounting for retention in roadside ditches, which delivered a result of roughly 0.60, i. e. 60 % of all PP emissions entering roadside ditches and subsequently storm drains are reaching surface waters via such bypasses.
With this introductory words said, the three factors used to define the case study scenarios are as follows: (i) transfer coefficients of 0.32, 0.46 and 0.60 corresponding to 54 %, 77 % and 100 % of 0.60 (TC0.32,TC0.46, TC0.60), (ii) discharge frequencies of one (T1) and six years (T6) (utilised for the calculation of the hydraulic radius/residence time and in turn PP retention) and (iii) the three described STREAM model flow direction scenarios, namely, topography (TOPO), tillage directions (TILL) and tillage directions combined with open furrows at all field borders (FURR).So all in all 18 scenarios (three transfer coefficients times two discharge frequencies times three flow direction data sets) were modelled and assessed.Storm drains at road embankments of about all asphalt roads were taken into account with the help of the storm drains model extension.A governmental reference routing dataset provided the necessary road data for this (Geoland.at, 2016).The period of the modelled scenarios ranged from the year 2008 to the year 2013.Water quality data of seven water quality gauges (additionally shown in Fig. 5) could be used to calculate mean annual PP river loads using the same methodology as of Hepp and Zessner (2019).
A detailed overview on all input data is provided by Zessner et al. (2016aZessner et al. ( , 2017)).Some other notable datasets nevertheless are (i) the already mentioned (geo-)database related to the IACS of the CAP of the EU contributing detailed information on the cultivated crops, the different factors of the USLE as well as the field borders (Hofer et al., 2014), (ii) a DEM with 10 × 10 m resolution, which served as the main input data for the STREAM model, (iii) non-agricultural land use based on the digital cadastre map and (iv) top soil characteristics derived from the digital soil map of Austria.The latter three datasets were supplied by the State Government of Upper Austria.Lastly, data on PP accumulation in top soil was extracted from Zessner et al. (2011Zessner et al. ( , 2016b) and Manning's roughness coefficients from Engman (1986).
Each of the nine T6 scenarios was then calibrated individually under a channel deposition rate of zero utilising the previously calculated mean annual PP river loads as targets and the overland deposition rate as the only calibration parameter.This more or less reflects long-term conditions including in-stream phosphorus stock depletion effects caused by major flood events (Zoboli et al., 2015) and can be considered a rather conservative approach for identifying critical source areas, as it generally leads to lower PP loads from cells farther away from surface waters due to higher overland deposition rates.The overland deposition rates of the nine T1 scenarios were subsequently adopted from the corresponding T6 scenarios thus simulating conditions as if no respective flood/transport event took place in the modelled period of six years.In order to achieve this, the discharge frequency related parameters of the hydraulic radius calculation are altered accordingly (Liu and De Smedt, 2004;Molnár and Ramírez, 1998).

Particulate PhozzyLogic Index (PPLI)
While providing interesting details, a map displaying the PP cell loads with, for example, 10 × 10 m resolution is not outright helpful to policy makers, who have to decide where to cost-effective implement a mitigation measure and where not.This information therefore has to be somehow transformed.For this purpose, we first calculate the total absolute PP contributions to surface waters for each field, i. e. the sum of all PP cell loads of each field with the help of zonal statistics and then apply a fuzzy membership function to those results.

Fuzzy membership function
A short introduction to fuzzy sets and fuzzy membership functions is given by Robinson (2003).Basically, fuzzy membership functions are functions turning data sets of arbitrary ranges and scales into fuzzy sets consisting solely of values between zero and one, where zero can be translated into "not possible" or not a member of a given set and one into "perfectly possible" or a member of a given set.All other values between zero and one represent a varying degree of possibility or membership in a given set.A value of, for instance, 0.7 thus stands for a higher degree of possibility or membership in a given set than a value of e. g. 0.4.At this point, it also has to be pointed out that possibility must not be confused with probability, even though they share the same value range.While probabilities are estimated based on data as well as certain assumptions and provide confidence intervals, possibility, especially due to its subjective choices of fuzzy membership functions, is less objective.Once one or more fuzzy membership functions have been chosen, it is perfectly objective nonetheless and can even be used to objectively compare expert judgements.This approach is particularly interesting and helpful in a sparse data environment as is the cost-effective implementation of mitigation measures against PP inputs into surface waters.
The fuzzy membership function we have chosen is the so called fuzzy large membership function and has the form where p 1 is the spread and p 2 is the midpoint.For the spread, we generally put in one, which makes the shape of the function somewhat similar to that of a logarithmic function, and for the midpoint, we put in a different value for each scenario.These values were calculated so that 80 % of the respective agricultural land received a value of less than 0.5 and 20 % a value of 0.5 or more (Pareto principle, cf.Sharpley et al., 2009).In this way, the 20 % of the total area of agricultural land belonging to the fields with the highest absolute PP inputs into surface waters received a value of 0.5 or more and are therefore what we consider critical source areas.Fig. 1 provides examples of fuzzy large membership functions with a midpoint of 4 kg yearly PP cell load per field and spreads ranging from one to five.With increasing spread, the function gradually transforms itself from a function with a shape similar to that of a logarithmic function to one with a shape similar to that of a sigmoid function.The black dashed lines shall help illustrate the meaning of the midpoint (possibility of 0.5).

Final PPLI creation
While the application of the described fuzzy large membership function to a single scenario already could be acknowledged a spatially distributed model result as well as fuzzy logic based PP index, it is only the first step, since fuzzy sets can further be overlaid in order to conclude.A simple technique for this is the so-called convex combination (Dubois and Prade, 1985;Kandel, 1986).In the process of convex combination, each involved fuzzy set can be assigned a weight, which makes it a weighted mean operator with weights adding up to one.Charnpratheep et al. (1997) also propose a modified version of this operator, which sets the overall result to zero in case one or more of the combined fuzzy sets are zero.
As no indication exists that one of our 18 scenarios is more possible than another as far as it comes to a single field and we neither favour one discharge frequency over the other, we used equal weights to overlay the fuzzy sets of all our scenarios with the unmodified convex combination operator in order to create the final PPLI.Using the modified operator does not contribute to the quality of the result in the absence of a knock-out criterion.

STREAM model flow directions
All in all, the flow direction of about 28 % (approx.440 000 cells) of the arable land were altered by enforcing tillage directions on topographic flow directions (TOPO vs. TILL).While comparing TOPO and FURR even results in a somewhat higher share of about 31 % (approx.500 000 cells), comparing TILL and FURR only shows a rather small share of about 8 % (approx.130 000 cells).This is reasonable, since the field borders merely comprise a small portion of the cells representing the arable land of the catchment.Additionally, these shares reveal that about 5 % of the cells were altered two times, first by the TILL and then by the FURR scenario.
Of the roughly 7000 arable fields nearly 6200 (approx.88 %) are affected by at least one cell with altered flow directions due to tillage directions and obviously all fields are affected by altered flow directions due to open furrows.
An example of the resulting flow directions of the three scenarios is given in Fig. 2.

Calibrated overland deposition rates
Calibration quality of the nine independently calibrated T6 scenarios is altogether very similar.They all exhibit a Nash-Sutcliffe efficiency (NSE) of about 0.95, a modified Nash-Sutcliffe efficiency (mNSE) of around 0.83, a percent bias (PBIAS) of almost zero and a ratio of the root mean square error to the standard deviation of measured data (RSR) of roughly 0.21 (Krause et al., 2005;Moriasi et al., 2007;Nash and Sutcliffe, 1970).Poesen (2018) considers scaling up sediment yields from plot to catchment scale one of the main challenges in  The calibrated overland deposition rates range from 0.91 to 1.46 × 10 −3 s −1 and decrease with lower transfer coefficients (1.16 to 1.46 × 10 −3 s −1 for 0.60 and 0.91 to 1.21 × 10 −3 s −1 for 0.32), which is reasonable, as lower inputs into surface waters from fields upstream of roads have to be compensated with higher inputs from fields upstream of surface waters.Generally, the deposition rates of the FURR scenario group are approximately 0.30 × 10 −3 s −1 lower than those of the other two STREAM model scenario groups.
These lower deposition rates of the FURR scenarios may appear counterintuitive.They are, however, a product of partially lower flow accumulations in the interior of the fields and higher flow accumulations along the field borders.
Apparently, the higher transport potential along the field borders affecting rather few cells does not fully compensate the lower transport potential affecting the comparatively many more cells of the field interiors.Thus, the process of calibration has to result in the overall lower deposition rates of the FURR scenario group in order to make up for the gap between the modelled and calculated PP river loads.

PP inputs into surface waters
Catchment-wide PP inputs into surface waters range from about 8.0 to 8.7 t yr −1 for the T1 scenario group and 15.9 to 16.2 t yr −1 for the T6 scenario group.Storm drains at road embankments account for approximately one quarter to almost one half of the PP inputs.The respective shares range from 24 to 29 % for the TC0.32 scenarios, 32 to 37 % for the TC0.46 scenarios and 39 to 44 % for the TC0.60 scenarios.This is in line with findings of Prasuhn (2011) whose ten year long field survey resulted in a share of 22 % of the total eroded soil.He also states that about half of the eroded soil was already retained within the borders of the source field though.This share therefore at least has to be doubled in order to account for the input into surface waters alone.Tab. 2 shows selected percentiles of the sums of the yearly PP cell loads per field of arable land for two selected scenarios.The distributions are as expected very skewed.In the case of the TC0.46-T1-TOPO scenario 70 % of the fields contribute less than 0.5 kg yr −1 PP to the total PP emissions into surface waters and even in the case of the TC0.46-T6-TOPO scenario the same percentile amounts to as little as nearly 1.1 kg yr −1 PP.The amounts of the higher percentiles are increasing rapidly from there on.

Identified critical source areas
The number of fields classified as a critical source area, i. e. fields with a possibility of 0.5 or more, ranges from 1151 to 1233 for all the modelled scenarios.14 538 fields, on the other hand, are not classified as a critical source area in any of the scenarios.This denotes an average of 1196 fields, which in relation to the total number of fields of 16 320 constitutes a share of about 7.3 %.The differences in the number of critical source areas can be regarded as negligible among the scenarios.
Furthermore, the median area of fields classified as a critical source area in at least one scenario is 2.94 ha with an interquartile range of 3.30 ha and the one of fields never classified as a critical source area is 0.74 ha with an interquartile range of 1.31 ha.These clear differences in the median field sizes and interquartile ranges are capable of explaining the relative small share of fields regarded as critical source areas compared to the approximately 20 % of the total area of agricultural land under consideration.
A major cause of this percentage mismatch is the choice of the sum of all PP cell loads per field as input to the fuzzy large membership function, which favours the selection of comparatively large fields.The purpose of this choice is the maximisation of the anticipated decrease in PP emissions into surface waters per field due to the implementation of suitable mitigation measures and, at the same time, the minimisation of the number of farmers involved.While other choices like the average of all PP cell loads per field may offer a better overall cost-effectiveness ratio, it may be beneficial to have fewer farmers involved, who have to be convinced to participate in a catchment-wide water quality protection programme.
The applied fuzzy large membership functions are designed to classify 20 % of the total area of agricultural land as critical source areas.Depending on the scenario, this share is responsible for 79 to 83 % of the total PP inputs into surface waters, which supports the 80:20 rule (Sharpley et al., 2009).Furthermore, apart from two outliers, all fields classified as a critical source area in at least one of the scenarios (1782 fields) are arable fields, hence the 20 % of the total agricultural land classified as a critical source area actually represent roughly 30 % of the total arable land.This nonetheless underlines the general principle at work.
Tab. 3 shows the number of times as well as the cumulative share in which the same field is classified as a critical source area.For the TC and T scenario groups the shares of fields classified as a critical source area in all of their respective scenarios lie despite their different group sizes around 50 %, whereas for the flow direction scenario groups (TOPO, TILL and FURR) the same shares account for somewhat more than 60 % of the fields.These differing shares can be explained by the fact that all flow directions and in turn flow accumulations are homogeneous within the flow direction scenario groups, but heterogeneous within the other scenario groups.So-called UpSet plots allow for a rather easy to interpret visualisation of multiple intersecting sets (Conway et al., 2017;Lex et al., 2014).Fig. 4 shows such a plot for a subset of the intersections of the critical source areas of all T6 scenarios.All nine scenarios agree on 816 (46 %) of the in total 1782 fields, which are classified as a critical source area at least once.Furthermore, there are 153 (96 plus 57; 9 %) critical source area fields where the FURR scenario group makes a clear difference.This once more emphasises the importance of knowing the relevant preferential flow pathways in a catchment under consideration.Another recognisable pattern is that the scenarios with lower transfer coefficients (especially TC0.32) differ slightly from those with higher transfer coefficients (TC0.46 and especially TC0.60).The T1 scenarios exhibit a similar pattern, albeit not as distinct as the one of the T6 scenarios.
All in all, the 18 scenarios perfectly agree on 720 of the 1782 fields classified as a critical source area in at least one of the scenarios.This is a share of about 40 %.The other 1062 fields or 60 % are diversely distributed among the various possible degrees of scenario intersections.Feasible policies now could be to only concentrate on the two fifth with perfect agreement or to implement mitigation measures on all of the 1782 critical source area fields.Including even more scenarios, however, could push these two numbers even further apart, particularly, since every of the assessed scenario groups comes with its own limitations.Table 3: Number of times as well as the cumulative share in which the same field is classified as a critical source area shown separately for each of the assessed scenario groups.Please note that the scenarios of the TC and flow direction scenario groups each occur only six times (the three scenarios of the respective other group times the two scenarios of the T scenario group), while the scenarios of the T scenario group each occur nine times (the three scenarios of the TC scenario group times the three scenarios of the flow direction scenario group).A limitation of the transfer coefficient scenario group is for sure the choice of a global transfer coefficient.This assumption may not hold, as areas with higher flow accumulations may also exhibit higher transfer rates.Another limitation is that the presence or absence as well as the exact location of individual storm drains remains uncertain (Hepp and Zessner, 2019).While the erosion part of PhosFate based on the (R)USLE reflects average conditions, its transport part has to be adjusted to a specific discharge frequency.With longer periods both parts eventually represent average conditions.Especially the T1 scenarios exhibit a mismatch in this regard.Increasing the discharge frequency (e. g. from six to one year) yet has a similar effect as increasing the overland deposition rate.Since the latter acts as calibration and therefore as catch-all parameter for, among others, the effects of catchment elements, which cannot be represented by the chosen spatial resolution of 10 × 10 m, but where retention can take place (e. g. unploughed strips between fields, hedges), the former can be used as a proxy for studying the effects of the presence or absence of such elements.The global nature of such an increase in discharge frequency once more poses a severe limitation, as it mainly influences long-distance transport.
Finally, the limitations of the flow direction scenarios are limited knowledge of the actual tillage directions, the surface roughness in and perpendicular to tillage directions and the exact locations of open furrows capable of concentrating flow.Some or even all of these shortcomings may be remedied by better available data or improved methods to derive them from existing data in the future.Tillage directions, for example, could be more accurately derived for L-shaped fields from field border data encompassing every single cultivated crop and not only the outer border of all crops in direct vicinity, which are cultivated by the same farmer.The latter, unfortunately, is the case with the field borders available for this study.

The final PPLI map
Fig. 5 presents the final PPLI map of the case study catchment.The fields coloured violet can be considered the catchment-wide most possible critical source areas.Most critical source area fields are located near the catchment boundary, where the slopes are comparably higher, and in the northwest of the catchment.This map, however, is the result of only a single option of overlaying multiple fuzzy sets.Further options, which should be tested in future research, are, for example, the application of different weights or other fuzzy overlay operators like the fuzzy gamma operator (Zimmermann and Zysno, 1980).
No matter how the different fuzzy sets are combined in the end, such a PPLI map allows each farmer to compare his or her fields to all the other fields within a certain catchment.This can be of great help when it comes to evaluating ifin the context of the entire catchment -a certain field substantially contributes to the overall PP emissions into surface waters or not.Furthermore, such a map provides a clear ranking for policy makers and could be used as a starting point for a state aid programme promoting the implementation of mitigation measures in a cost-effective way.For this purpose, the height of subsidies could even be linked to the degree of possibility of a field being a critical source area or not.
Yet we consider it important that the final decision on the kind and exact location of the implemented mitigation measure stays with the farmer, as he or she has the best local knowledge.For example, in many cases only a small portion of vegetated buffer strips receives the majority of overland flow (Djodjic and Villa, 2015;White and Arnold, 2009), hence, the concept of critical source areas can even be translated to a single field.Allowing farmers to choose from a variety of mitigation measures and combine them in a smart way may be capable of significantly improve their implementation quality and in turn effectiveness.

Conclusions
The presented algorithm for the allocation of PP emissions entering surface waters to their cells of origin is a major step in the direction of a model based and, in future, potentially automatic identification of critical source areas on catchment scale.As a straightforward single flow algorithm using PP transports as weights while taking care that mass balances are not violated, it has the potential to work also together with multi-flow algorithms.Switching to a multi-flow algorithm may be necessary in case it is desirable to increase spatial resolution even further.
Further important puzzle pieces for an automation of the identification of critical source areas are zonal statistics in conjunction with fuzzy logic.Utilising zonal sums, fuzzy large membership functions with a spread of one and the convex combination operator with equal weights all in all worked well for the creation of the final PPLI, however, future research should be carried out in order to assess if there are more advantageous choices for each of these steps and to develop improved or new methods for an automatic derivation of tillage directions as well as especially the specific positions of open furrows.In addition, depending on catchment characteristics, need may arise to specify different scenarios or even apply different models (e. g. physically based).As long as there is a way to calculate per-field statistics from the model results, they are eligible for the PPLI.

1.
Initialising an intermediate cell load either as the transported amount of PP entering a surface water (riparian cells -first cycle) or as the apportioned cell transfer in the last step of the previous cycle (all other cells -subsequent cycles).2. Initialising an intermediate cell transfer by setting it equal to the previously initialised intermediate cell load.3. Updating the intermediate cell load by multiplying it with the ratio of the local and the local plus the sum of the inflowing PP transports.4. Setting the final cell load utilising the local net PP emission as upper limit (minimum of the previously updated intermediate cell load and local net PP emission).5. Calculating the cell load carry-over as the maximum of the intermediate cell transfer minus the final cell load and zero.6. Setting the final cell transfer to the previously calculated cell load carryover.

Figure 1 :
Figure 1: Examples of fuzzy large membership functions with a midpoint of 4 kg yearly PP cell load per field and spreads ranging from one (dark blue curve) to five (yellow-green curve).With increasing spread, the function gradually transforms itself from a function with a shape similar to that of a logarithmic function to one with a shape similar to that of a sigmoid function.The black dashed lines shall help illustrate the meaning of the midpoint (possibility of 0.5).

Figure 2 :
Figure 2: Example of the resulting flow directions of the TOPO (top left), TILL (bottom left)and FURR scenario (bottom right).Especially the FURR scenario leads to partially more extreme flow accumulations, as it concentrates flow at field borders and limits spillovers to neighbouring fields to cells with a single possible downstream path.

Figure 3 :
Figure 3: Individual yearly PP cell loads (left) and transfers (right) of a small area within the case study catchment for the TC0.46-T6-TOPO scenario.

Figure 4 :
Figure 4: UpSet plot for a subset of the intersections of the critical source areas of all T6 scenarios.While the set sizes of the nine sets are presented in the bottom left and their intersection sizes in the top right area, the bottom right area provides information about the sets involved (black dots connected with a black line or single black dots) in each of the above displayed intersections.
Final PPLI map of the case study catchment: fields with an overall possibility of 0.5 or more, which can considered the most possible critical source areas, are coloured violet.A darker tone thereby indicates a higher possibility.

Table 1 :
Additional properties of the case study catchment.

Table 2 :
Selected percentiles of the sums of the PP cell loads per field of arable land in kg yr −1 for two selected scenarios.