Coastal margins and backshores represent a major sink for marine debris: insights from a continental-scale analysis

Marine debris represents a major threat for the environment. Plastic production is increasing exponentially and causing an unprecedented growth of plastic pollution entering the marine environment. Hence, a thorough assessment of debris accumulation areas is required to address the longstanding question about where is all the missing plastic. Most research on marine debris sinks to date has focused on oceanic gyres, the water column, seabeds and wildlife. Relatively little has focused on the potential of coastal areas as debris sinks. To address this knowledge gap, the spatial distribution of debris from the waterline to the backshore was modelled from a continental-scale dataset of coastal debris distribution from 635 surveys across Australia. Results showed that the distribution of debris is significantly correlated with oceanic and atmospheric processes (i.e. onshore Stokes drift and wind), and coastal usage for recreational activities (i.e. regional population and distance to the nearest road). Debris density and size increased from the waterline to the backshore, indicating that the backshore area represents an important debris sink, especially for larger sized items.


Introduction
The increased abundance and persistence of marine debris in the environment represents an unprecedented and urgent ecological, economic and societal issue. As the magnitude of the problem has grown over the last decades, so has the engagement of the scientific community and the awareness by the general public. Considering the negative impacts marine debris have on the environment, a global consensus has developed on the need for a more circular economy and improved waste management solutions (UNEP 2016;Ellen MacArthur Foundation 2017, Velis et al 2017.
Among the different types of material that comprise marine debris, plastic is regarded as the most harmful to the environment. Consequences span from biodiversity losses (Derraik 2002, Gall andThompson 2015) and negative impacts on the development and reproduction of many wildlife species (Oehlmann et al 2009, Foley et al 2018, Bucci et al 2020 to potential threats to human health (Thompson et al 2009, Smith et al 2018, World Health Organisation 2019 and reduction in tourism revenue (Jang et al 2014, Leggett et al 2014. Global plastic production is doubling approximately every 11 years (Plastics Europe 2013), and reached 359 million tonnes in 2018 (Plastics Europe 2019). This means that we will produce as much plastic between now and 2030 as we have since production began in 1950. The durability of plastic items, regarded as one of the major strengths of the material, has become a tremendous downside given its mismanagement and loss into the environment. This is particularly relevant considering the fast pace at which global plastic production is increasing. At the same time, plastic pollution in the marine environment is growing rapidly (Lebreton et al 2018, Wilcox et al 2019. It was estimated that 4.8-12.7 million metric tonnes of mismanaged plastic waste entered the ocean in 2010 alone, with the Asia-Pacific identified as the most significant leakage hotspot (Jambeck et al 2015, Lebreton et al 2017. The African continent, undergoing unprecedented population growth and urbanisation, has the potential to turn into another major source of land-based waste entering the ocean if robust waste management policies and solutions are not implemented and enforced in the coming years (Jambeck et al 2018, Lebreton andAndrady 2019).
Different model estimates show that the amount of plastic particles afloat on the surface of the ocean is more than 5 trillion pieces (Eriksen et al 2014), ranging between 15 and 51 trillion pieces weighing between 96 and 236 thousand metric tonnes (Van Sebille et al 2015). As the latter model estimate accounts for approximately only 1% of the amount of mismanaged plastic waste predicted to enter the ocean in 2010 (sensu Jambeck et al 2015) and the rate of plastic accumulating in the ocean gyres is increasing in time (Lebreton et al 2018, Wilcox et al 2019 the question arises about where the remaining 'missing' plastic waste is located (Thompson et al 2004).
While plastic has been reported on beaches around the world for more than four decades (Scott 1972, Cundell 1974, most of the research on 'missing' marine debris has focussed on the oceanic gyres (Lebreton et al 2018), the water column (Dai et al 2018), the bottom of the ocean (Woodall et al 2014, Chiba et al 2018, or seabirds and wildlife (Boerger et al 2010, Schuyler et al 2014, Roman et al 2016. Recently however, Brennan et al (2018) argued that coastal environments need to be included as possible sinks for anthropogenic debris. Estimates of yearly riverine inputs of plastic waste span from coastal populations are between 1.15 and 2.41 million metric tonnes (Lebreton et al 2017). These inputs from surface waters moving through urban areas result in higher deposition of plastic debris along the coast near outflows, suggesting that much of the material may be deposited along the coast near its input to the ocean (Willis et al 2017), explaining the mismatch between estimates of loss from land and those of standing stock and rates of increase offshore.
To investigate this mismatch between estimates of loss from the land boundary and changes in the load of plastic in the sea, we focused on understanding the processes that determine the spatial distribution of debris in the coastal environment, from the waterline to the backshore area. Following the analysis performed by Hardesty et al (2017), we used a continental scale dataset on debris loads on the Australian coastline to investigate (1) the amount of debris on the coast, (2) its spatial distribution with respect to the marine-terrestrial boundary, (3) the role of debris size in determining its position between the waterline and the terrestrial environment, and (4) the evidence for oceanic and atmospheric onshore transport of marine debris. These questions were assessed by developing two separate survival regression models, which also allowed us to draw conclusions on the role of the global coastline as a sink for mismanaged plastic lost into the ocean.

Field surveys
Field data were collected along the coast of Australia between September 2011 and June 2016 (figure 1), following the methodology presented in Hardesty et al (2017). In brief, surveys were conducted in all Australian coastal states and territories at approximately 100 km intervals. Due to access limitations, the northern part of the country had significantly fewer sampling sites than the rest of the Australian coastline. After a random selection of the first survey site in the northeastern part of Australia, sites were selected approximately 100 km apart, with a minimum of three and maximum of 6 transects surveyed at each location (see detail below). At each sampling site the access point location was recorded via GPS receiver, together with the date, time of day, weather conditions (including wind direction and a categorical measurement of wind speed), name of the observer, and number of visible people. When the information was available, the amount of time since the last beach cleanup was recorded for the site. The characteristics recorded for each transect included: Start and end location, length, position of the dominant debris line, start and end time of the survey, shape and aspect of the coastline, beach gradient, substrate type and colour, and backshore type. At each site, the first transect was located at least 50 m from the access point and each subsequent transect was located 50 m from the preceding transect. For each transect, a tape measure was laid down extending from the water's edge to two metres into the backshore vegetation (or seawall). A two-metre belt was considered for the counts of debris (one metre per side, with one observer per side). The pieces of debris counted during the field surveys represent the standing stock of debris, either stranded or directly deposited, hence representing equilibrium conditions. Items had to be detected from standing height in order to be recorded. The total length of the transect, from the waterline to two metres into the backshore was divided in ten equal length sections (figure 2). In this study the term waterline indicates the line dividing dry sand from the coastal area affected by active wave motion, while the term backshore indicates the area located behind the coastal site where no evidence of sustained marine influence was visible, with backshore types ranging from terrestrial vegetation (i.e. trees, shrubs, grass, etc) to cliffs, dunes and anthropogenic seawalls (see SI table 1 is available online at stacks.iop.org/ ERL/15/074037/mmedia). In the case of vertical cliffs and seawalls, the transects ended at the location of the cliff or seawall. For each of the ten sections, the distance to the first piece of debris was recorded, together with its size (on a categorical scale; see figures 3(B), (C), material and colour. The latter two characteristics were also recorded for every piece of debris found along each transect.

Physical and anthropogenic forcing data
In addition to the data collected in situ, for each site we obtained GIS data of population counts in a circular area from each survey site within radii of 5, 20, and 50 km; and distance from the access point to the closest road (Australian Bureau of Statistics 2016). To assess the magnitude of physical forcing, measurements were obtained for wind speed, significant wave height and mean wavelength for the entire year before the sampling date from three different forecast models. Atmospheric data were obtained from the National Centre for Environmental Prediction (NCEP) Climate Forecast System Version 2 (CFSv2) 6 hourly Products (Saha et al 2011a(Saha et al , 2011b and Selected Hourly Time-Series Products (Saha et al 2011a(Saha et al , 2011b. Ocean wave data were obtained from CAWCR Wave Hindcast (Durrant et al 2013). Significant wave height and mean wavelength were used to calculate the shoreward Stokes drift velocity at the 20 m isobath (Bever et al 2011). Stokes drift was relevant for this study as it represents the net drift velocity of a particle in the direction of wave propagation (Van den Bremer and Breivik 2017).

Statistical analysis and modelling
All statistical analyses were performed in R i3.5.0 environment (R Core Team 2018). The goal was to understand the factors affecting the measured distances to the first debris item in each of the 10 sections of each coastal transect, noting that the distance to the first piece of debris in a transect section is inversely proportional to debris density in that section (Hilborn and Mangel 1997). Because transects were different lengths, positions were converted to a relative scale where location was measured as the proportion (in each of the 10 sections) from the seaward to the landward bound of the transect (as shown in figure 2). Because debris could be present or absent within each section, survival analysis ('survival' package; Therneau 2015) was utilised to model the distance covered from the beginning of the section until the first item was found (Crawley 2012). When no debris was found in a section (N=0), the observation was treated as right censored, as per standard procedure in survival analysis. In this case, right censoring indicates that the distance at which debris appeared is greater than the length of the section, but otherwise unknown. To include this in the model, a dummy variable status was created and each section was assigned a value of either 0 or 1 reflecting the absence or presence of the first debris item, respectively. As the instantaneous risk of finding a piece of debris was considered to be constant along the section, the data was fit with an exponential distribution (Crawley 2012).
Two statistical models were developed to perform the analysis. One comprised the entire dataset, and hence is referred to as total, while the other included only the sections where debris was found, and is referred to as positives. The main difference between the two, other than sample size, was the inclusion of censored observations in the total model. These pose some limitations on the variables that could be incorporated in the model, as only variables that were directly observed for all sections, regardless of their status, were allowed. Consequently, as the role of debris size could not be investigated in the total model, a reduced version of the dataset was analysed following the same distribution used for the total model to assess whether size played a significant role in the spatial distribution of debris when present in the coastal environment.
To evaluate the effects of site-specific characteristics and external forcing, the covariates included in the models were divided into three main groups: (1) observation process; which included aspect, backshore, gradient, debris size category (for the positives (C) Size category of the first items found in sections 1, 5 and 10 (size category is indicated by the colour, consistent with 2B). Letters above the bars correspond to distributions that are (letters different) or are not significantly different (letter shared) at the p<0.05 level based on a Komogolov-Smirnov test of the shape of the distributions. model), coastal shape, state, substrate colour, substrate type, total debris count in the transect, and weather; (2) physical forcing; which included mean onshore wind velocity and mean onshore Stokes drift velocity integrated over timescales varying between one day and half a year; and (3) anthropogenic forcing which included scaled population within a circle with radius of 5, 20, and 50 km, and distance from the access point to the nearest road. The observations group included covariates that have an impact on debris retention and deposition (but not position, e.g. coastline shape) based on past analyses , Brennan et al 2018 and that may affect debris detection (e.g. substrate colour). The covariates belonging to the physical forcing group were included in the model to test for the impact of wind and wave motion on the deposition of debris from the ocean and its subsequent transport on land, while those belonging to the anthropogenic forcing group were included to assess the contribution of potential terrestrial inputs.
The best overall total and positives models were obtained following analogous procedures. For each step, the best model was selected among all the possible combinations of covariates included using the second-order Akaike Information Criterion (AICc), as calculated from the dredge function ('MuMIN' package; Bartoń 2018). The AICc value of a model, analogously to the AIC, is a relative indication of the fit of the model to the data. Ranking models based on AICc score, as done by the dredge function, allowed for robust selection of the best model, because the AICc score takes into account the parsimony of each model (Burnham andAnderson 2002, Symonds andMoussalli 2011). Initially, the predictive capacity of the covariates in the observation group was tested. Once the best observation model was identified, two new models were implemented: One including the best observation model and the physical forcing covariates, and the other including the best observation model and the anthropogenic forcing covariates. The final model was identified by combining the covariates included in the best versions of these two forced models, and assessing the AICc of all possible combinations.

Results
Analyses were performed on 635 transects (6350 sections), distributed across 188 sites located around Australia (figure 1). The minimum transect length was 2 m, while the maximum was 1860 m (median=30 m, mode=20 m). Less than 15% of the transects were characterised by cliffs and anthropogenic seawalls as backshore types, while dunes accounted for 34% of the transects. Three locations (1.6%), 144 transects (23%), and 4996 sections (79%) were characterised by the complete absence of debris. The presence of debris increased steadily from the waterline to the backshore, with section 1 representing only 4% of those sections where debris was observed and section 10 representing 18% of the sections where debris was present ( figure 3(A)). The most abundant materials observed along transects were plastics (56%; of which 73% was hard plastic, 16% was soft plastic, and 11% was plastic string), glass (17%), foam (10%), rubber (5%), metal (3%) and paper (3%).
Almost half of the first items found in each section belonged to the smallest classes: class 1 (0.2-1.2 cm) accounted for 20% of the total and class 2 (1.2-2.4 cm) for 25% ( figure 3(B)). No items with any dimension larger than 36.3 cm were found as first items. The size distribution of first items found in sections 1, 5 and 10 highlights the differences in debris characteristics from the waterline (section 1) and in the backshore area (section 10; figure 3(C)). Size distributions in these three sections were compared using the Komogolov-Smirnov test, which confirmed a significant difference (p=0.03) in size distribution between section 10 and the other two, at locations 1 and 5.
According to AICc values of the total model (table 1), among the covariates that belong to the observation group, those included in the best observation model were aspect, backshore type, section number, state, substrate colour, substrate type, total debris and weather (AICc=5844.9). Of the different timescales assessed for mean onshore wind velocity and mean onshore Stokes drift velocity, the most significant were one week for the wind and one month for Stokes drift (AICc=5751.9). On the other hand, the anthropogenic forcing covariates that gave the best outcome when added to the best observation model were population within a 50 km radius, distance to the nearest road, and the interaction term between these two (AICc=5761.6). The best full model consisted of all the covariates included in the best observation plus physical forcing model and population within a 50 km radius, the distance to the nearest road, and the interaction between the section number and the population within a 50 km radius terms (AICc=5650.5, table 1).
Mean onshore Stokes drift and all the anthropogenic forcing covariates encompassed in the best full model were significantly related (p<0.05) to the relative distance at which the first piece of debris was located (table 2; full version in table 1S). The observation covariates that had a stronger relationship with the position of the first piece of debris were section number and total debris in the transect. According to the direction of the coefficient and its relative size with respect to the median, section number and mean onshore Stokes drift are the strongest predictors in terms of effect size, and distances to the first item decreases with both, suggesting the backshore during periods of heavy wave activity has the shortest distances, and thus the highest densities.
When reducing the dataset to the positives model (table 3), the best observation model (AICc=631.7) included debris size, section number, total debris and the interaction term between debris size and section number. This was a better model than those in which physical forcing was included. In contrast, the incorporation of covariates belonging to the anthropogenic forcing group allowed for an improvement (i.e. decrease) of the AICc score when the population within a 50 km radius term was included (AICc=627.1). As no covariates within the physical forcing group were significant, and the inclusion of the interaction term between section number and population within a 50 km radius gave a higher (i.e. worse) AICc value, the best full model corresponded to the best observation plus the anthropogenic forcing model. The significant terms (p<0.05) in the best positives full model (table 4) were total debris, population within 50 km, and the interaction term between size and section number. As total debris and population within 50 km increased, the relative distance from the beginning of the section to the first piece of debris decreased. Interestingly, the coefficient of the interaction term between size and section number was positive, indicative of an increase in relative distance when both size and section number increased. Model predictions of the relative distance at which debris was found in each section based on its size showed that variations in relative distance between the lower and upper sections of the beach were more accentuated for larger debris (figure 4).

Discussion
Initially posed by Thompson et al (2004), a longstanding question has remained about where the 'missing' plastic in the ocean is located. This uncertainty has persisted for more than a decade in spite of the recent increase in studies addressing plastic in the ocean. To date, researchers have favoured the notion that the plastic is located at the bottom of the ocean Table 1. AICc scores for the total model. For each group, all the models within 2 AICc units from the best model are reported, as they fall within the 95% confidence interval around the best model, denoted by the lowest AICc. Factors included were as follows: A=aspect, B=backshore, G=gradient, Sn=section number, S=state, Sc=substrate colour, Ss=substrate type, Td=total debris, W=weather, S1Dy=mean onshore Stokes drift over one day, S1Mo=mean onshore Stokes drift over one month, W1Wk=mean onshore wind speed over one week, W10Dy=mean onshore wind speed over ten days, WHyr=mean onshore wind speed over half a year, P50=population living within a 50 km from the site, Rd=distance from the access point to the closest road.

Model
Factors AICc  . However, few have questioned whether indeed the plastic is actually missing. This work and results from Brennan et al (2018) suggest that there is substantial amount of debris deposited on the coastline. Data and modelled results obtained from this study showed that the density and size of debris located in the coastal environment increases from the waterline to the backshore. This provides strong evidence that backshore areas should be considered as marine debris sinks. Moreover, the size distribution of first items found in the sections ( figure 3(B)) showed that stranded debris detected at the study sites is significantly larger than debris located in the upper ocean layer in samples offshore of the Australian continent, where median size is 2.8 mm (Reisser et al 2013; see Morét-Ferguson et al 2010 for data on the North Atlantic). Even if the median size of the items found in the upper ocean layer is much smaller than that found in the lowest and middle sections of the beach, the significant difference in the distribution of item sizes observed in sections 1 and 5 and section 10 may be explained by the dominance of ocean-sourced, smaller items in the lower sections of the beach, which are the most influenced by wave action.
Aspect, backshore, substrate colour and type were found to be associated with the relative distance at which debris was located. Their effect was in general agreement with a previous analysis that assessed the quantity of debris in the coastal environment (Hardesty et al 2017). Indeed, factors that were associated with an increase in total debris along the transect in the model proposed by Hardesty et al (2017) were negatively correlated to relative distance in this study, as one would expect. As discussed by Hardesty et al (2017), these variables largely describe factors that affect either the probability of detection of plastic during surveys (e.g. substrate colour) or overall rates of retention of debris that is deposited (Brennan et al 2018). For the purposes of understanding the spatial distribution of items between the water boundary and the terrestrial boundary, however, they are important to account for, but are not the key factors of interest.
Both onshore wind and Stokes drift were found to have a significant relationship with debris density. While the two processes act in the same direction (i.e. shoreward), the effect of Stokes drift was substantially larger than that of onshore wind (table 1). The different timescales on which onshore wind and onshore Stokes drift act, one week and one month respectively, potentially separate the processes of debris stranding (Stokes drift) and redistribution (wind). Not only onshore wind would redistribute debris stranded in the lower sections of the beach, but stronger wave action would allow ocean-sourced debris to strand closer to the backshore during stronger storms. Onshore wind exposure and backshore type were previously identified as critical factors in determining the resuspension rate of debris (Brennan et al 2018). Therefore, the higher density of debris found in the upper part of the coast may be the result of Stokes drift and onshore wind transport towards the backshore and a low departure rate from this area of accumulation once material is deposited there. This is illustrated in figure 2, where one can envision it would be difficult for debris items to escape from the vegetation at the terrestrial boundary once deposited there, in line with findings by Brennan et al (2018).
The number of people living within 50 km from a survey site was negatively correlated to the distance at which debris was found in the sections along a transect, reflecting the positive correlation between regional population and debris density in the coastal environment (sensu Hardesty et al 2017). We found this effect to be particularly strong for sections closer to the waterline. Noting that the total load of debris was already included as a separate term in the two survival models, this suggests that ocean-sourced and directly deposited debris are found at different locations at a coastal site. The higher debris density closer to the waterline was accentuated when the regional Table 3. AICc scores for the positives models. For each group, all the models within 2 AICc units from the best model are reported, as they fall within the 95% confidence interval around the best model, denoted by the lowest AICc. Factors follow from table 2. In addition, Sz=debris size; W1Dy=mean onshore wind speed over one day; W1Mo=mean onshore wind speed over one month; S10Dy=mean onshore Stokes drift over ten days.

Model
Factors AICc population was larger, highlighting direct deposition of debris closer to the shore, potentially driven by people spending more time there. In contrast, oceansourced, stranded debris appears to be mostly in the upper sections and in the backshore (figure 2). We suggest that direct deposition by beach users is not the only source of debris in the lower coastal sections, but it superimposed on the process of debris stranding driven by Stokes drift and onshore wind, before debris is picked up by onshore wind and pushed toward the backshore.
Not surprisingly, more accessible beaches were found to have more debris in each section, pointing again to the underlying mechanism of increased numbers of people on beaches affecting debris densities and distribution via littering, as opposed to ocean transport from nearby population centres. The effect of regional population on debris density in the coastal environment highlights the necessity of waste management regulations and policies, especially in areas where the population is larger. Identifying this as a significant factor driving debris at coastal sites leaves room for optimism considering that regional governments and populations have the opportunity to make local changes, and suggesting that local changes can lead to significant reductions in local debris density (Willis et al 2018).
The distance at which the first item was detected in each section was found to be associated with the size of the item itself (figure 4), indicating the predominance of smaller debris in the lower sections and of larger debris in the upper sections of each transect. This may be explained by the different characteristics of distinct parts of the coastline. Section 1, the part of the transect closest to the water, together with the lower sections, is affected more by wave motion at the waterline. The constant physical interaction between wind, waves, and substrate likely causes debris to be broken into smaller pieces, which are more responsive to wave motion and may move inland and shoreward frequently with ongoing wave action. In contrast, given its cross-sectional area and lower relative density, larger debris may experience greater forcing, driving it further onshore than smaller items. Indeed, more and larger debris was located further from the waterline (figure 3(C), section 10; figure 4). Different to the rest  of the transect, the last two metres of section 10 consist of backshore, which has a significant effect in decreasing departure rates of debris from the transects (Brennan et al 2018). Thus, the increase in distance at which larger debris is found in section 10 may reflect backshore accumulation, especially of larger items. While it is possible that large debris in the backshore vegetation are observed due to ease of detection compared to smaller debris in vegetated areas and that the presence of vegetation might obscure smaller items, the thorough sampling approach taken, consisting of two experienced observers per transect recording debris from standing height whilst moving at a consistent rate, means this is unlikely to be the reason for observed difference in item size. Consistently, small material is still detected in the backshore, although at lower rates than in the middle of the transect, which is generally near the main debris line on the beach.

Conclusion
While numerous studies have addressed plastic in the ocean and quantified debris in coastal environments, previous work had not asked where along the coastline debris density is highest. The development of two survival regression models based on coastal surveys that quantified debris from the waterline into the coastal zone, allowed for the identification of factors that drive the overall distribution of debris from the waterline to the backshore and enabled us to challenge the longstanding question of where is all the missing plastic. Mean onshore Stokes drift over one month was found to be the dominant oceanic process, while mean onshore wind over one week was the dominant atmospheric process associated with coastal debris density. The effect of onshore Stokes drift resulted to be much larger than that of onshore wind in determining the relative distance of debris, and the effects of wind and Stokes drift varied on debris of different sizes. Evidence was found to suggest that larger regional population was associated with higher debris density in each section. This effect progressively weakened from the lowest to the highest section, suggesting that direct deposition of debris by beach users in the lower sections of the beach was prevalent, and that larger regional populations and more access were associated with more debris near the waterline.
Lastly, the amount and size of debris located in the coastal environment increased from the waterline to the backshore, suggesting that backshore area, particularly when vegetated, should be regarded as a major sink for marine debris, in particular for larger items. Therefore, future studies assessing the quantity of 'missing' plastics (Thompson et al 2004) should include estimates of debris located in backshore areas.