Modelling Beach Litter Accumulation on Mediterranean Coastal Landscapes: An Integrative Framework Using Species Distribution Models

Beach litter accumulation patterns are influenced by biotic and abiotic factors, as well as by the distribution of anthropogenic sources. Although the importance of comprehensive approaches to deal with anthropogenic litter pollution is acknowledged, integrated studies including geomorphologic, biotic, and anthropic factors in relation to beach debris accumulation are still needed. In this perspective, Species Distribution Models (SDMs) might represent an appropriate tool to predict litter accumulation probability in relation to environmental conditions. In this context, we explored the applicability of a SDM–type modelling approach (a Litter Distribution Model; LDM) to map litter accumulation in coastal sand dunes. Starting from 180 litter sampling plots combined with fine–resolution variables, we calibrated LDMs from litter items classified either by their material type or origin. We also mapped litter accumulation hotspots. LDMs achieved fair-to-good predictive performance, with LDMs for litter classified by material type performing significantly better than models for litter classified by origin. Accumulation hotspots were mostly localized along the beach, by beach accesses, and at river mouths. In light of the promising results achieved by LDMs in this study, we conclude that this tool can be successfully applied within a coastal litter management context.


Introduction
The presence of macro litter on sandy beaches and dunes can cause health problems [1,2], economic losses [3,4], and damage to natural biodiversity [5][6][7][8]. In light of that, investigating litter distribution and accumulation is a key issue in coastal zone sustainable management [9]. Currently, the amount of waste washed up on coasts is a matter of great concern in Mediterranean countries, and waste mismanagement has been considered one of the most important environmental problems affecting the coastal courtiers of Europe [6,10,11]. The European Commission adopted the Marine Strategy Framework Directive (MSFD, 2008/56/EC), a normative instrument aimed to achieve a 'good environmental status' of European marine and coastal environments [12]. Under the MSFD, Member States are committed to undertake specific actions related to research, monitoring, and surveillance programs to keep "the properties and quantities of marine litter below a threshold above which it could cause harm to the coastal and marine environment" [12].
Marine litter can be categorized according to its source (e.g., households, agriculture, industrial activities, recreational activities, fisheries, shipping [10]) and to its composition (e.g., plastic, glass, paper, aluminum, polystyrene [13]). Overall, the distribution of litter types classified by source is supposed to be mainly related to the presence of anthropogenic structures and human activities [14,15]. On the contrary, the distribution of litter items grouped by material types likely depends on their physical characteristics, e.g., buoyancy and resistance to seawater, as well as their susceptibility to be blown by the wind [9,16].
Several studies looked into the geographical patterns of beach litter [9] and suggested that they are influenced by abiotic and biotic drivers, such as, e.g., vegetation structure [17,18], dune morphology [19], wind, wave action, tides [20], and the density of debris materials [18,21,22]. On the other hand, other studies suggested that marine litter accumulation is mostly influenced by anthropogenic sources [23], identifying polluted rivers [24][25][26], densely populated urban areas [17,27], seashore mass tourism [25,28], fisheries, shipping, and thriving aquaculture activities [29][30][31] to promote litter accumulation on the coasts. Although the role of these factors in driving litter occurrence is well known, and the importance of a comprehensive approach for dealing with anthropogenic litter pollution has been recognized [12], only a few studies have investigated the simultaneous effect of these drivers on litter accumulation patterns on the coasts. Moreover, studies integrating geomorphologic, biotic, and anthropic factors are needed to better understand and predict accumulation patterns of beach litter [15,21,25]. In this regard, most of the studies investigating beach litter transportation and accumulation patterns focused on very complex numerical modelling approaches (e.g., [32][33][34][35]), which somehow hamper their applicability within a context of coastal zone management aimed at fulfilling MSFD recommendations.
In the fields of ecology and conservation biology, a class of modelling approaches called "Species Distribution Models" (SDMs) [36] designed to quantify the occurrence probability of a species in a given area has been successfully used by scientists and environmental managers (e.g., [37][38][39]). More recently, SDMs were successfully applied in non-strictly ecological contexts. In fact, such a technique has been used to predict the occurrence probability of roadkill events for wildlife (e.g., [40,41]) or the susceptibility of a given region to landslides [42]. Considering these examples, SDMs might analogously represent an effective tool for predicting litter accumulation probability in relation to given environmental conditions. Being based on statistical relationships between the occurrence of a given event (dependent variable; here, litter fragments) and geographical predictors (explanatory variables), SDMs could be suitable for both assessing litter-environment relationships and for predicting their spatial distribution. Furthermore, quantifying the relationship between observed spatial patterns of litter occurrence and different geographical variables (e.g., geomorphology, vegetation, anthropogenic factors, etc.) can enhance our knowledge about the process of litter accumulation.
In this perspective, we proposed a study aiming to explore the applicability of a SDM-type modelling approach (a Litter Distribution Model; LDM) to map beach litter accumulation in coastal dunes accounting for the simultaneous effect of geomorphologic features, biotic conditions, and anthropogenic pressure. In doing that, we assume litter accumulation is not homogeneous across the coastal landscape but it occurs wherever the combined influence of dune morphology, vegetation pattern, and anthropic pressure is favorable to the litter occurrence. Within this framework, the specific objectives of our study were to (i) calibrate LDMs for litter items classified by source and material type; (ii) compare LDM predictive accuracy for these two main classification schemes; (iii) explore which factors are associated with beach litter accumulation hotspots; and (iv) map such accumulation hotspots.

Study Area
The study was conducted along 88 km of the Adriatic coastline in central Italy (Abruzzo and Molise regions; Figure 1), considering six protected areas where Mediterranean coastal dune vegetation is well represented [43]. The analyzed coasts include long, continuous sandy beaches and recent (Holocene) dunes usually occupying a narrow strip along the sea-shore, where the deposition of sediments and organogenic materials is mostly due to wave motion and wind [44,45]. Under natural conditions, these dunes are characterized by a complex sequence of ecosystems that follow the sea-inland ecogeomorphological gradient [46,47], ranging from herbaceous annual plant communities on the strandline lower zone of the beach, crossing through perennial herbaceous vegetation on embryonic and shifting dunes, to patchy Mediterranean scrubs on the inland stabilized dunes, and to Pinus sp. woods in the foredunes [48,49]. As for most of the Mediterranean coasts, several tracts of the study area are seriously impacted by human pressure that promotes both land-based [50,51] and ocean-based litter production, movement, and deposition [17,19]. Since the mechanical cleaning of beaches affects the integrity of natural dune eco-morphology [43], beach litter removal in protected areas is done with low frequency [16]. In light of that, protected areas offer an ideal scenario to investigate which environmental, geomorphological, and anthropogenic factors may drive beach litter accumulation.

Study Area
The study was conducted along 88 km of the Adriatic coastline in central Italy (Abruzzo and Molise regions; Figure 1), considering six protected areas where Mediterranean coastal dune vegetation is well represented [43]. The analyzed coasts include long, continuous sandy beaches and recent (Holocene) dunes usually occupying a narrow strip along the sea-shore, where the deposition of sediments and organogenic materials is mostly due to wave motion and wind [44,45]. Under natural conditions, these dunes are characterized by a complex sequence of ecosystems that follow the sea-inland eco-geomorphological gradient [46,47], ranging from herbaceous annual plant communities on the strandline lower zone of the beach, crossing through perennial herbaceous vegetation on embryonic and shifting dunes, to patchy Mediterranean scrubs on the inland stabilized dunes, and to Pinus sp. woods in the foredunes [48,49]. As for most of the Mediterranean coasts, several tracts of the study area are seriously impacted by human pressure that promotes both land-based [50,51] and ocean-based litter production, movement, and deposition [17,19]. Since the mechanical cleaning of beaches affects the integrity of natural dune eco-morphology [43], beach litter removal in protected areas is done with low frequency [16]. In light of that, protected areas offer an ideal scenario to investigate which environmental, geomorphological, and anthropogenic factors may drive beach litter accumulation.

Beach Litter Occurrence Data
Beach litter sampling was conducted through a stratified random protocol based on 2 × 2 m survey units [18,27]. Specifically, we randomly placed 180 plots considering, as sampling strata, the EU habitat types, occurring in the study area (RanVegDunes database [48]; Natura 2000 habitat maps, ftp://ftp.dpn.minambiente.it/Natura2000/TrasmissionECE_ 2013/schede_mappe). This sampling protocol and plot size are particularly effective for describing the highly heterogeneous eco-geomorphology and the fine-scale environmental gradients that characterize coastal dunes [47,52]. Sampling was carried out in the vegetative period (April-May 2018), making sure no mechanical cleaning had been carried out in the study area for several months [30]. All the plots were visited only once by a single team of researchers. Litter items > 2.5 cm (i.e., macro litter [53]) were collected, visually inspected, and recorded following the OSPAR protocol (Convention for the Protection of the Marine Environment of the North-East Atlantic) [54]. The position of each sampling unit was georeferenced using a GPS in order to ensure future monitoring campaigns [30]. Subsequently, we classified litter items according to two macro categories, i.e., material type and origin, following the OSPAR protocol and its successive integrations proposed for the Mediterranean area [55] (Table 1). We calibrated a different LDM for each material type and origin category. Frequently, occurrence data used in species distribution modelling come from online databases and/or citizen science initiatives. Accordingly, such data may be spatially biased and need to be appropriately filtered (e.g., [41]). Since our modelling exercise relies on litter data gathered through a statistically robust sampling design, we avoided implementing any ex-ante data filtering, even though we found it appropriate to test for spatial autocorrelation in LDM residuals [39] in all cases (see below).

Geographical Covariates
We hypothesized that the occurrence of accumulation areas can be affected by a set of geomorphologic, biotic, and anthropogenic predictors. Starting from airborne LIDAR imagery (acquired in 2008 by the Italian Ministry of Environment, within the PST-A mission), we derived two geomorphological predictors, i.e., curvature and slope, along with a canopy height model (CHM). Indeed, dune morphology (concavity, convexity, height) and vegetation (forests, shrublands or herbaceous plant communities) may likely play a role in influencing litter movement as well as in capturing litter items transported by wind [7,16], thus regulating spatial patterns of accumulation areas. In addition, we considered a number of predictors related to human activity, such as the Euclidean distance from artificial surfaces [56] (derived from the imperviousness products of the Pan-European High Resolution Layers of the Copernicus service; https://land.copernicus.eu/pan-european/highresolution-layers/imperviousness), from harbors, and direct access to the beach [17,27]. We also included the Euclidean distance from river mouths, hypothesizing a contribution to accumulation areas by litter items transported by rivers [25,27]. The position of harbors, accesses, and river mouths was manually digitalized and their Euclidean distance from each survey unit was calculated. All the predictors were rasterized at a spatial resolution of 2 m (i.e., the native spatial resolution of LIDAR imagery) and checked for their multicollinearity by posing a variance inflation factor ≤ 3 [57]. All the variable preparation and selection steps were carried out within the R environment [58].

Litter Distribution Models
We calibrated the LDMs by using the maximum entropy modelling algorithm implemented in MAXENT version 3.3.3k [59]. This algorithm compares covariate values at site localities to a sample of background points to create a map of probability of event occurrence ranging from 0 (i.e., no accumulation probability, in this case) to 1 (highest accumulation probability). Since MAXENT predictions are sensitive to initial modelling settings [60], we tested different MAXENT implementations through the ENMeval R package [58] to find the settings that optimize the trade-off between goodness-of-fit and overfitting [61]. In particular, we tested regularization values between 0.5 and 4, with 0.5 steps. Furthermore, we included the following feature classes: linear, linear + quadratic, hinge, linear + quadratic + hinge, linear + quadratic + hinge + product, and linear + quadratic + hinge + product + threshold [41,61]. Among the resulting 48 combinations, we chose the one reporting the lowest Akaike Information Criterion (AICc) [62]. To calibrate LDMs, a set of 10.000 background points was randomly placed within the six investigated protected areas to describe their environmental and geomorphological characteristics and to represent pseudo-absences. We evaluated LDMs through a spatial block cross-validation scheme relying on the inherent data partition structure due to the protected areas network. Specifically, we calibrated LDMs with data from five out of six protected areas, leaving the held-out data for evaluation purposes. We repeated the procedure by holding out in turn the data from each protected area. The block cross-validation scheme proved able to assess model transferability, i.e., the ability to extrapolate predictions into new areas [63] and to penalize models based on biologically meaningless predictors [64]. The predictive performance of the models was assessed by measuring the area under the receiver operating characteristic curve (AUC [65]), the difference between calibration and evaluation AUCs (AUCdiff [66]), and the true skill statistic (TSS [67]). AUC values range from 0 (models with no predictive ability), to 1 (perfect predictions [68]), while TSS values range from −1 (no predictive ability) to 1 (perfect prediction [67]). We predicted LDMs over the six protected areas included in the analysis and applied the "ExDet" approach [69] to assess the effect of model extrapolation on values of predictor variables lying outside the calibration range (i.e., negative values of the D metric indicate extrapolation). LDM projections were binarized according to three threshold approaches (i.e., "maximize TSS", "minimum training presence", and "10th percentile" [70]) to account for the effect of using different binarization schemes [71]. For each litter class, binary maps obtained under the three thresholding schemes were stacked and summed [72] to obtain spatially explicit predictions of litter accumulation hotspots.

Hotspot Analysis
The spatial association between litter accumulation hotspots, i.e., areas where multiple litter classes are predicted to occur, and the geographical covariates considered in LDMs was explored by calibrating a Random Forest (RF) algorithm, where the pixel-by-pixel count of each litter source/material obtained by summing binary maps was used as the response variable. We also included the different thresholding schemes as additional categorical covariate to account for the different hotspot maps generated by the three thresholds considered. As RF parameters, we set (i) 1000 decision trees, (ii) number of variables randomly selected at each node equal to two (i.e., the square root of the number of geographical covariates), and (iii) the Gini index as the split rule [73,74]. RF goodness-of-fit was assessed through the out-of-bag R 2 (i.e., the mean prediction error of each RF training sample x i , using only the trees that did not include x i in their bootstrap sample [75]). To ease evaluating the association strength between litter accumulation hotspots and geographical covariates, we followed the approach proposed by Igras and Biecek [76] and fitted spline functions through the RF marginal response for each covariate. Then, such splines were used as surrogate covariates in generalized linear models (GLMs) against litter pixel-by-pixel count (i.e., as done for RF models). The statistical significance of each surrogate covariate coefficient quantified in these GLMs allows us to evaluate if the association between litter hotspots and geographical covariates as predicted by RF is strongly supported by the data [76].

Litter Distribution Models
LDMs achieved fair-to-good predictive performance (sensu Swets [68] and Landis and Koch [77]), showing AUC values ranging from 0.765 (SD = 0.203) for "packaging" source to 0.874 (SD = 0.130) for "aluminum" material type; AUCdiff values scored from 0.060 (SD = 0.178) for the "food and beverage" origin to 0.035 (SD = 0.150) for the "other" origin, and TSS values ranged from 0.551 (SD = 0.101) for "plastic" material type to 0.844 (SD = 0.085) for "aluminum" material type (see also Table S1). Although LDMs for both litter classification schemes yielded good predictive performances, models for litter classified by material type reported significantly higher TSS values than models for litter classified by origin (Figure 2; Wilcoxon test W = 583, p < 0.05). Negligible spatial autocorrelation was found in LDM residuals (Appendix A), and only ca. 11% (SD = 0.031%) of the study area reported a marginal extrapolation effect (ExDet D value = −0.034, SD = 0.013; Figures S1-S10). Most of the litter classes based on both material type and origin classification exhibited similar responses to a specific pool of geographical covariates that reported a variable contribution >10% in almost each LDM. Among the most recurrent responses, litter accumulation probability increased toward low CHM values (e.g., containers, fishing, mixed materials, packaging), close to direct accesses to the beach (e.g., aluminum, plastic) and near to river mouths (e.g., polystyrene, glass), and far from artificial surfaces (e.g., food; Figures S11-S20). Most of the litter classes based on both material type and origin classification exhibited similar responses to a specific pool of geographical covariates that reported a variable contribution >10% in almost each LDM. Among the most recurrent responses, litter accumulation probability increased toward low CHM values (e.g., containers, fishing, mixed materials, packaging), close to direct accesses to the beach (e.g., aluminum, plastic) and near to river mouths (e.g., polystyrene, glass), and far from artificial surfaces (e.g., food; Figures S11-S20).

Hotspot Analysis
Overall, accumulation hotspots hosting the highest number of litter classes (i.e., five) were mostly localized along the beach and embryonic dunes (Figure 3), as well as near direct accesses to the beach and to river mouths ( Figure 3E,F; see also Figures S21-S25). Hotspot distribution was not homogeneous within each protected area. For instance, protected areas A and B were largely covered by continuous hotspots potentially hosting all the material types ( Figure 3A,B). On the contrary, protected areas E and F include small and isolated hotspots ( Figure 3E,F). Considering litter classes by material type, most of the study area was covered by hotspots hosting four litter classes (ca. 700 ha), while hotspots with five litter classes covered ca. 300 ha (Figure 4). Considering litter classes by origin, hotspots with five classes were the widest, covering more than 1000 ha in the study area, whereas the other hotspots occupied around 250 ha ( Figure 4).   RF models achieved an excellent goodness-of-fit for both material type and origin hotspots, reporting out-of-bag R 2 values of 0.974 and 0.965, respectively. For both litter classification criteria, RF models highlighted CHM and distance from river mouths and direct accesses to the beach as the most important covariates driving accumulation hotspot occurrence ( Figure 5 and Figure S26). Specifically, the count of litter classes increased close to vegetation with a low canopy height (e.g., shrubs) and near to river mouths and direct accesses ( Figure 5). GLM results indicated all the surrogate covariates fitted through the RF marginal response were significantly associated with the count of the litter classes (Tables S2 and S3), suggesting strong statistical support of the RF results.  RF models achieved an excellent goodness-of-fit for both material type and origin hotspots, reporting out-of-bag R 2 values of 0.974 and 0.965, respectively. For both litter classification criteria, RF models highlighted CHM and distance from river mouths and direct accesses to the beach as the most important covariates driving accumulation hotspot occurrence ( Figure 5 and Figure S26). Specifically, the count of litter classes increased close to vegetation with a low canopy height (e.g., shrubs) and near to river

Discussion
In the present study, we showed how a modelling approach such as SDM can be effectively used to map beach litter accumulation patterns in a given area. In fact, LDMs achieved good accuracy levels in predicting beach litter occurrence along the Central Adriatic coastal dune ecosystems, explicitly integrating the influence of geomorphologic features, biotic conditions, and anthropogenic pressure in a single analytical framework. Furthermore, we showed that classifying litter items by their material types yielded more accurate LDMs than models calibrated from litter grouped by origins. Interestingly, our results highlighted that the influence of geographical variables on litter accumulation patterns varies across the coastal landscape. Specifically, debris occurrence probability increased toward dune vegetation with low canopy height values, close to direct accesses to the beach and to river mouths, and far from artificial surfaces, thus appearing jointly driven by context-specific biotic conditions, geomorphologic features, and anthropogenic factors. A similar effect was observed for plant canopy height and distance from river mouths and direct accesses in shaping litter accumulation hotspots.

Alternative Litter Classification Schemes Generate Different Predictive Accuracies
We found LDMs from litter classified by material type (e.g., plastic, polystyrene, glass, paper, aluminum, etc.) achieved good predictive performance according to both AUC [68] and TSS [77], suggesting their occurrence patterns are closely associated with the geographical covariates included in the models. This evidence might likely be related to the specific buoyancy and weight of each material type, which influence the transport dynamics mediated by wave and wind, resulting in similar accumulation patterns. Along the Mediterranean coasts, waves and tides usually drop off beach litter along the drifting line, embryonic, and mobile dunes [17,27]. These act as a sort of "source area" from where litter can be blown farther inland [7] until perennial herbaceous vegetation of shifting dunes and shrubs of fixed dunes stop and entangle the litter [26]. On the other hand, heavy materials tend to remain in place and be trapped in/under sand of the shifting dunes [78,79]. Geographical covariates, e.g., geomorphological features and plant canopy height, which we included in LDMs, were able to approximate such dynamics, therefore adequately explaining geographical variation in litter occurrence patterns and leading the models to perform well. While LDMs from litter classes by origin achieved overall acceptable performance (fair-to-good according to Swets [68] and Landis and Koch [77]), their accuracy was significantly lower than that for models for litter classes by material type. While the influence of the source on the litter accumulation pattern at the coarse scale is well documented in the literature [23,30], classifying litter items according to their origin puts a less direct focus on their buoyancy and weight, which likely represent the characteristics that mostly interact with geographical covariates. As a result, the modelled relation between the occurrence of litter classes by origin and such covariates was weaker. That said, we cannot entirely exclude that the lower accuracy of litter origin LDMs might rather depend on an improper split of the litter classes within the origin criterion.

Geographical Factors Affecting Litter Accumulation Hotspots
In the Mediterranean dunes system analyzed in this study, biotic elements such as vegetation height showed a predominant role in shaping litter accumulation patterns, confirming similar evidence documented by previous research [7,56]. Actually, litter accumulation patterns seem to follow the sea-inland gradient that characterizes the Mediterranean coastal dunes [17]. In particular, litter occurrence was strongly related to low-medium values of canopy height (e.g., containers, fishing, mixed materials, packaging), indicating a higher accumulation across the foredune vegetation zonation. Most of the heavy litter remains in the aphytic zone and on the strand line [80], as well as a small part of light litter that is trapped by organic waste as woody posts [7]. Moving toward inner sectors occupied by perennial herbaceous vegetation on embryonic and shifting dunes, litter items tend to be partly entangled by grasses and tufts, whose root structure is particularly able to fix and stabilize sediments [81,82], thus leading to litter accumulation and burial. Light litter items (e.g., polystyrene) can be further blown up toward fixed dunes, where most of them are curtailed by the pioneer scrubs of Mediterranean maquis [27]. Therefore, our results highlighted as a novel finding that not only the seashore and embryonic dunes but also foredunes, with perennial herbaceous and woody vegetation, are hotspots of litter accumulation across low-lying coastal areas. The perennial vegetation works as a "flying plastic litter sink" and a biological barrier that greatly reduces the accumulation of these items in the tall pine and evergreen oak woodlands growing in the adjacent inner sectors. Such an "uncommon" pattern is likely related to the dune characteristics in the study area. In fact, closed seas as the Adriatic are more protected from prevailing westerlies than open sea basins. Therefore, coastal dunes form "low-lying" systems [45] that are particularly prone to a sparse litter accumulation pattern, able to invade multiple dune sectors. Unfortunately, the litter removal in dense grasses and pioneer maquis is too hard to do and therefore polystyrene and plastic bags trapped under the branches become persistent and cumulative over time [16]. Further, the "flying plastic litter sink" vegetation is included in habitats of European concern (Directive 92/43/CEE) and is inhabited by several animal species included on the Red List such as the Kentish plover (Charadrius alexandrinus [83]) and the Hermann's tortoise (Testudo hermanni [84]).
As similarly observed along oceanic coasts [85], on Mediterranean seashores, the river discharge into the sea contributes to beach litter accumulation by providing a variety of litter items coming from inland alluvial plans. In fact, the prevalent occurrence of litter hotspots near river mouths evidenced by our results is coherent with well-established literature reporting litter accumulation to be favored near coastal rivers [19], which act as pathways of terrestrial mismanaged waste to the sea [25,85]. Litter items discharged by rivers can then be reshuffled and deposited on the seashore by sea currents and tides, therefore remaining mostly close to river mouths [34].
The direct accesses to beaches emerged as another important determinant of beach litter accumulation hotspots, as well as of several single litter items (e.g., aluminum, food, and plastic; see Figures S11, S14 and S18). Such rather intuitive evidence is supported by several studies showing how small roads and trails increase the potential number of visitors frequenting the seashore, which in turn increase litter droppings [86]. It is widely documented how recreational activities represent one of the main factors behind beach litter accumulation [25,30,87]. Their effects appear even worse in protected natural areas due to an overall poor ecological awareness of the tourists. Moreover, in such areas, it often necessary to walk long distances before reaching the closest cans to deposit garbage [29,88].
Although not resulting among the most relevant covariates, we reported an interesting effect of the distance from urban areas. While several studies in the Mediterranean basin reported a higher marine litter abundance close to populated and industrialized areas [89], as well as a clear relationship between growing populations in coastal cities and increasing coastal debris accumulation [90], our results showed an inverse relationship (see e.g., Figures S11-S20). Interestingly, more recent papers found similar evidence [25,87,91], interpreting the lower litter amount close to the urban areas as due to an efficient litter collection by the municipal litter agencies. Such evidence might be valid in our case, although we did not have detailed information about the cleaning effort in the study area to confirm this.
Contrary to other similar studies [92], we found no relevant effect played by geomorphological features (i.e., slope and curvature) in shaping litter accumulation patterns. Actually, this seems a reasonable outcome given the characteristics of the costal dune system analyzed in our study. Since such a system is mostly shaped by low-lying dunes [45], it is plausible that its mild geomorphologic features did not have a relevant effect in driving accumulation patterns. That said, we could not exclude a priori a role of geomorphology also in our peculiar context. That is the reason why we have included, and suggest to consider, dune morphology among the covariates in LDMs.

Management Implications
The integration of geomorphologic features, biotic elements, and anthropogenic pressure into LDMs allowed us to better investigate litter occurrence patterns and yielded accurate spatially-explicit predictions of litter accumulation hotspots along the analyzed coast. Incorporating ecological modelling tools into a coastal zone management context offered a new possibility of mapping areas with high pollution hazard, which can support site-specific management actions. In fact, a better understanding of the effect played by the individual factors behind litter accumulation, along with the production of debris accumulation maps, represent essential steps for developing effective strategies coherently to those committed to by the MSFD. From such a perspective, the strong geographical variability in litter accumulation patterns in the protected area network analyzed in this study supports the approach of prioritizing the implementation of tailor-made measures to control, monitor, and prevent the formation of litter hotspots, in order to achieve good coastal environmental status. For instance, in the analyzed low-lying coastal areas which are particularly prone to storm surges [45] that wash up important amounts of litter reaching embryonic dunes, periodic manual cleaning of the driftline is advisable. The observed hotspots of items blown by the wind to fixed dunes and trapped on perennial vegetation evidenced a paucity of such periodical cleaning campaigns. In such cases in which waste collection is occasionally possible, cleaning work must be preferentially done in litter hotspot areas.

Conclusions
The modelling framework implemented in this study represents a new, promising tool in the context of coastal zone management, able to explicitly integrate several factors notoriously involved in litter accumulation patterns. In keeping with an ever wider tendency of exporting SDM toward different fields, we explored if and how this statistical tool, which is highly popular among ecologists, might perform in predicting litter accumulation occurrence. We cannot exclude that adding some specific predictors (e.g., waves, tides or currents) involved in litter transport, especially in the water, might have improved the obtained predictions. Unfortunately, we did not have mapped information on such parameters at a sufficiently high detail along the study coastline to be included in our models. Therefore, we assumed a preponderant effect of dune morphology, vegetation pattern, and anthropic pressure on the modelled litter accumulation patterns. While we are aware of this basic assumption, we found it acceptable especially given our specific focus on terrestrial litter accumulation (i.e., dunes and river mouths). That said, it would surely be worthwhile to adapt our approach in future exercises focusing on marine accumulation. Despite such limitations, our study showed that, when fed with highly accurate litter occurrence data combined with geographical covariates at extremely high resolution (e.g., LIDAR products), LDMs are able to predict litter accumulation with good predictive performance. In addition, this approach allowed exploring the differential role in litter accumulation played by factors such as geomorphology, vegetation, and human infrastructures, in terms of their magnitude and shape. Different from other, very complex approaches relying on numerical modelling, the outputs of LDMs can widely support preventive and effective measures in litter hotspots to avoid further litter accumulation in coastal EU habitats, where it is detrimental for wild plants and animals and its removal is extremely difficult. In light of the promising results achieved by LDMs as described in this study, we conclude that this tool can be successfully extended to different coastal zone management applications. Notwithstanding, we stress the importance of replicating our approach in different coastal area contexts, in order to further test its reliability.
Supplementary Materials: The following are available online at https://www.mdpi.com/2073-4 45X/10/1/54/s1, Figure S1: Map of ExDet metric values showing possible extrapolation of LDM predictions for aluminum; yellow to red tones indicate low to high extrapolation, while blue tones refer to no extrapolation. Figure S2: Map of ExDet metric values showing possible extrapolation of LDM predictions for containers; yellow to red tones indicate low to high extrapolation, while blue tones refer to no extrapolation. Figure S3: Map of ExDet metric values showing possible extrapolation of LDM predictions for fishing and boating; yellow to red tones indicate low to high extrapolation, while blue tones refer to no extrapolation. Figure S4: Map of ExDet metric values showing possible extrapolation of LDM predictions for food and beverage; yellow to red tones indicate low to high extrapolation, while blue tones refer to no extrapolation. Figure S5: Map of ExDet metric values showing possible extrapolation of LDM predictions for mixed materials; yellow to red tones indicate low to high extrapolation, while blue tones refer to no extrapolation. Figure S6: Map of ExDet metric values showing possible extrapolation of LDM predictions for other materials; yellow to red tones indicate low to high extrapolation, while blue tones refer to no extrapolation. Figure S7: Map of ExDet metric values showing possible extrapolation of LDM predictions for packaging; yellow to red tones indicate low to high extrapolation, while blue tones refer to no extrapolation. Figure S8: Map of ExDet metric values showing possible extrapolation of LDM predictions for plastic; yellow to red tones indicate low to high extrapolation, while blue tones refer to no extrapolation. Figure S9: Map of ExDet metric values showing possible extrapolation of LDM predictions for polystyrene; yellow to red tones indicate low to high extrapolation, while blue tones refer to no extrapolation. Figure S10: Map of ExDet metric values showing possible extrapolation of LDM predictions for glass; yellow to red tones indicate low to high extrapolation, while blue tones refer to no extrapolation. Figure S11: Variable importance (bar plot) and response curves describing the shape of the relationship between aluminum accumulation probability (y axis) and each explanatory variable (x axis); red curves refer to the mean response, while blue ribbons indicate ±1 SD around the mean. Figure S12: Variable importance (bar plot) and response curves describing the shape of the relationship between container accumulation probability (y axis) and each explanatory variable (x axis); red curves refer to the mean response, while blue ribbons indicate ±1 SD around the mean. Figure S13: Variable importance (bar plot) and response curves describing the shape of the relationship between fishing and boating accumulation probability (y axis) and each explanatory variable (x axis); red curves refer to the mean response, while blue ribbons indicate ±1 SD around the mean. Figure S14: Variable importance (bar plot) and response curves describing the shape of the relationship between food and beverage accumulation probability (y axis) and each explanatory variable (x axis); red curves refer to the mean response, while blue ribbons indicate ±1 SD around the mean. Figure S15: Variable importance (bar plot) and response curves describing the shape of the relationship between mixed material accumulation probability (y axis) and each explanatory variable (x axis); red curves refer to the mean response, while blue ribbons indicate ±1 SD around the mean. Figure S16: Variable importance (bar plot) and response curves describing the shape of the relationship between other material accumulation probability (y axis) and each explanatory variable (x axis); red curves refer to the mean response, while blue ribbons indicate ±1 SD around the mean. Figure S17: Variable importance (bar plot) and response curves describing the shape of the relationship between packaging accumulation probability (y axis) and each explanatory variable (x axis); red curves refer to the mean response, while blue ribbons indicate ±1 SD around the mean. Figure S18: Variable importance (bar plot) and response curves describing the shape of the relationship between plastic accumulation probability (y axis) and each explanatory variable (x axis); red curves refer to the mean response, while blue ribbons indicate ±1 SD around the mean. Figure S19: Variable importance (bar plot) and response curves describing the shape of the relationship between polystyrene accumulation probability (y axis) and each explanatory variable (x axis); red curves refer to the mean response, while blue ribbons indicate ±1 SD around the mean. Figure S20: Variable importance (bar plot) and response curves describing the shape of the relationship between glass accumulation probability (y axis) and each explanatory variable (x axis); red curves refer to the mean response, while blue ribbons indicate ±1 SD around the mean. Figure S21: Distribution map of the accumulation litter material-type hotspots identified by LDM projections binarized by the "minimum training presence" threshold. Numbers and colors indicate the cumulative number of litter classes (e.g., plastic, polystyrene, paper, glass, aluminum and mixed materials). Figure S22: Distribution map of the accumulation litter material-type hotspots identified by LDM projections binarized by the "10th percentile" threshold. Numbers and colors indicate the cumulative number of litter classes (e.g., plastic, polystyrene, paper, glass, aluminum and mixed materials). Figure S23: Distribution map of the accumulation litter origin hotspots identified by LDM projections binarized by the "maximize TSS" threshold. Numbers and colors indicate the cumulative number of litter classes (e.g., containers, fishing and boating, food and beverage, packaging, other materials). Figure S24: Distribution map of the accumulation litter origin hotspots identified by LDM projections binarized by the "minimum training presence" threshold. Numbers and colors indicate the cumulative number of litter classes (e.g., containers, fishing and boating, food and beverage, packaging, other materials). Figure S25: Distribution map of the accumulation litter origin hotspots identified by LDM projections binarized by the "10th percentile" threshold. Numbers and colors indicate the cumulative number of litter classes (e.g., containers, fishing and boating, food and beverage, packaging, other materials). Figure S26: Bar plot depicts the relative contribution of the different covariates in driving the accumulation hotspots of litter origin. Line plots report the shape of the relationship between the count of litter classes in hotspots and the three most influential covariates. Red lines refer to the marginal response as predicted by RF models, while blue lines indicate the spline curves fitted through such responses and used as surrogate covariates in GLMs. Table S1. Evaluation metrics reported by LDMs for each litter element. The table includes mean ± standard deviation values referring to cross-validation replicates. Table S2. Coefficients of generalized linear models including surrogate geographical covariates fitted through the random forest marginal response. The RF model was calibrated for litter material-type classification. Table S3. Coefficients of generalized linear models including surrogate geographical covariates fitted through the random forest marginal response. The RF model was calibrated for litter origin classification. Data Availability Statement: Data available on request due to restrictions eg privacy or ethical. The data presented in this study are available on request from the corresponding author. The data are not publicly available due to ongoing research activities.

Acknowledgments:
This study was carried out with the partial support of the Interreg Italia-Croazia CASCADE (CoAStal and marine waters integrated monitoring systems for ecosystems proteCtion AnD management-Project ID 10255941) and PON-AIM (Programma Operativo Nazionale ricerca e innovazione 2014-2020; ID AIM1897595-2). The authors acknowledge the editor and the anonymous reviewers for their comments that helped us to improve the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Checking for Significant Spatial Autocorrelation in LDM Residuals
Appendix A. 1

. Method
We calculated the models' residuals as 1-predicted probability of the presence for each litter item record, also including the pseudo-absences. Moran's I was calculated considering multiple distances between points, ranging between a minimum distance with each point connected only to its nearest neighbour to a maximum with all points connected. Significance of Moran's I was calculated using a randomization test with 999 Monte Carlo permutations. All procedures were repeated 10 times, each time randomly choosing 1000 pseudo-absences.

Appendix A.2. Results and Conclusions
All of the litter items, classified either by origin or by material type, exhibited a significant but weak correlation. In particular, Container litter reported the highest correlation, showing a mean Moran's I value = 0.147 ± 0.501, which was statistically significant in ca. 90% of the replicates. By contrast, Food and beverage litter scored a mean Moran's I value = 0.006 ± 0.200, which was statistically significant in ca. 75% of the replicates. Averaged among all the analysed litter items, we reported a mean Moran's I value = 0.009 ± 0.259, which was statistically significant in ca. 70% of the replicates. Such low values allow considering as negligible the degree of correlation of LDM residuals.