## 3 Results
### Benefits of NbS to Mitigate Soil Erosion Risk  
We seek to find the NbS that decreases the component of the shear stress responsible for soil erosion \(\tau_{bs}\) below the critical threshold \(\tau_{c}\). The main objective is to impede that \(\tau_{bs}>\tau_{c}\) as occurs actually in most of the olive-grove area (simulation CN89-A described in Section 3.1 and Figure 6c). Otherwise, the significant shear stresses provoke unsustainable soil losses and the formation of deep gullies in most of the catchment (e.g., Figure 7a). If we achieve the goal \(\tau_{bs}<\tau_{c}\), we stop the actual erosion processes and reduce the catchment sediment budget. Furthermore, we prevent the silting of the Guadalen Reservoir, where the Manillas stream drains at a distance of 4 km from the basin outlet. Sediments exposed to erosion over the Manillas basin are mostly coarse silt, with a uniform grain size in the range of \([0.031,0.0625]\) mm. The critical shear stress for the onset of sediment motion was determined from Table 7.1 in Julien [77], which yields \(\tau_{c}=0.083\) N\(\cdot\)m\({}^{-2}\) in water at 20 \({}^{\circ}\)C. For evaluating the local shear stress responsible for erosion \(\tau_{bs}\), we used the classical shear-stress partitioning that decomposes the effective boundary shear stress \(|\mathbf{\tau}_{b}|=\tau_{bf}+\tau_{bs}\) into the form drag \(\tau_{bf}\) and the skin friction \(\tau_{bs}\) (see details in SS 2.8.1 by [PERSON] [54]). A simple and compact way to visualise the benefits of NbS in the whole watershed is to analyse the changes in the probability of the flow depth (Figure 10a), velocity magnitude (Figure 10b), and shear stress (Figure 10c) regarding the current situation. To compute the probabilities, we accounted only for the regions of interest: the olive groves and the drainage network (blue and red areas in Figure 4a, respectively). Hence, we excluded the non-erodible urban areas and civil infrastructure from the analysis. In the absence of an NbS (solid black lines in Figure 10), the most probable values of the flow depth, velocity, and shear stress vary in the narrow intervals \([1.6,2.4]\) mm, \([0.3,0.4]\) m\(\cdot\)s\({}^{-1}\), and \([3,5]\) N\(\cdot\)m\({}^{-2}\) for all the infiltration rates (simulations CN100-A, CN89-A, and CN64-A). Indeed, the three curves in black nearly overlap. The wetland (NbS D in purple) and the reservoir pond (NbS E in red) do not provoke a substantial modification of the flow depth, which is consistent with the observed observed (Figure 10). Figure 10: Histograms of the simulated (**a**) flow depth \(h\), (**b**) magnitude of the velocity vector \(|\mathbf{u}|\), and (**c**) the skin friction \(\tau_{\text{bs}}\) for the fifteen scenarios explained in Section 3.3.  
the probability curves that are similar to the actual ones (scenario A in black). A significant part of the shear-stress curves for the simulations type A, D, and E lie out of the target region \(\tau_{bs}<\tau_{c}\). Table 2 shows the cumulative probability for the catchment area where erosion occurs due to \(\tau_{bs}>\tau_{c}\). Independently of the infiltration rate, erosion develops in more than 86.4% in scenarios A, D, and E. The values in Table 2 show that the infiltration does not play a role in controlling erosion because of the same reasons that did not substantially attenuate the streamflows (Figure 9a). The steep slopes of the watershed and the high precipitation intensity favour runoff formation to the detriment of infiltration processes. The interventions on the agricultural area (NbS C, green lines in Figure 10c), including the channels of the drainage network (NbS B, blue lines in Figure 10c), changed the shear-stress distribution by a great extent. The distribution switched to a multimodal distribution. The major mode decreased from \([3,5]\) N\(\cdot\)m\({}^{-2}\) (present-day) to \([2.4,6.5]\times 10^{-4}\) (NbS B) and \([5.1,9.2]\times 10^{-4}\) (NbS C) when the mean curve number varied in the range of \([64,100]\).