## Abbreviations  
1. Use of [PERSON]'s et al. [13] methods to optimize the W2 WQMP: non-metric multi-dimensional scaling (NMDS), principal component analysis (PCA), Kriging, Moran's Index and leave - one - out cross-validation. First two multivariate analysis were applied: NMDS: add to the coordinate of each sampling station the distance to the outflow of the watershed and a PCA to associate to each station one or two factorials created from the concentrations of respectively PT, CF, and TSS for the complete data series (2009-2016). The amount of data (n) and the relation between the variables confirm the prerequisite of PCA, that of multi-normality and the linearity of the relations. Outliers are a source for distortion of these tests and were removed. Considering the missing values for PT, CF and TSS we chose not to reduce the data series, but we chose to select reconstruction of the missing data with the package MissMDA with R--for further information consult [26; 27]. Second, we proposed several Kriging models ([30] with the objective to identify the sampling stations which are contributing the most to predict the spatial model of the watershed and the distribution of these stations. The quality of the models was judged on the comparisons of the standard error of the residues, the values of the adjusted R\({}^{2}\), the significance of the model and the criteria of Akaike (Akaike information criterion--an estimator of relative quality of statistical models). We tested six models, and for each one and we chose to retain the following model:  
\[\text{Model M6}:\quad Log[FC]=Log[TSS]\times Log[TP]\times\text{{Dist}}\]  
where FC = Average Fecal coliform concentration, TSS = Average Total suspended solids concentrations, TP = Average Total phosphorous concentrations and Dist = distance in linear kilometers between each sampling station and the outflow of the watershed (W2). For each model, we used the spatial distribution of the residues of the leave-one-out cross-validation to visualize the increase in the variance of the residues of the model when they are removed. These stations are then considered as important, as their removal harms the model. On the other hand, stations for which the retrieval does not harm the model (cause little or no variance) are considered as those for which the information is redundant and which are suggested for removal [8]. Finally, we also calculated for each pair of stations the Moran Index and the associated Z score. After several tests, we retained the method called Minimum Spanning Tree, which proposes a structure were the sum of the distances between stations is minimal. Some relations were modified in order to obtain a representative relationship between the stations and the hydrographic network. The Z scores based on the Moran Index calculated on this structure provides an indicator of the redundancy of information between two or more neighboring stations. For each station, a score Z \(\leq\) -1,96 translates to a strong dissimilarity between a station and its neighbouring stations. A score of Z \(\geq\) 1,96 indicates a strong redundancy of the values of one station and its neighbouring stations [29]. For the optimization objective, verify whether sampling frequency was sufficient, we used the general linear regression for each of the variables for the entire sampling period (2009-2016) for the seven stations listed in Table 1 (method tested by [16]). In order to obtain a sampling effort of 50 %, we created two subsets of data for each station by choosing one date out of two from the original data set. To obtain a sampling effort of 33.33%, we chose one date out of three from the original data sets to obtain three data sets. The subsets were modelized and presented as linear regressions lines superposed on the model including all the available data.