Modelled transport of benthic marine microplastic pollution in the Nazaré Canyon

With knowledge of typical hydrodynamic behavior of waste plastic material, models predicting the dispersal of benthic plastics from land sources within the ocean are possible. Here we investigated the hydrodynamic behavior (density, settling velocity and resuspension characteristics) of non-buoyant preproduction plastic pellets in the laboratory. From these results we used the MOHID modelling system to predict what would be the likely transport and deposition pathways of such material in the Nazaré Canyon (Portugal) during the spring/summer months of 2009 and the autumn/winter months of 2011. Model outputs indicated that non-buoyant plastic pellets would likely be transported up and down canyon as a function of tidal forces, with only a minor net down canyon movement resulting from tidal action. The model indicated that transport down canyon was likely greater during the autumn/winter, primarily as a result of occasional mass transport events related to storm activity and internal wave action. Transport rates within the canyon were not predicted to be regular throughout the canyon system, with stretches of the upper canyon acting more as locations of pellet deposition than conduits of pellet transport. Topography and the depths of internal wave action are hypothesized to contribute to this lack of homogeneity in predicted transport.


Introduction
Marine microplastic pollution in the forms of preproduction pellets, fragments, filaments, films and foams originates from direct spillage and breakdown of plastic debris (Moore et al., 2011) and synthetic materials with densities ranging from ∼ 0.9-1.4g mL −1 (Morét-Ferguson et al., 2010;Andrady, 2011).Plastics of higher density (HD) than seawater are concentrated in marine and fluvial benthic environments (Galgani et al., 2000;Claessens et al., 2011;Costa et al., 2011;Mordecai et al., 2011).Data on sub-surface microplastic abundance and distribution within the marine environment is limited in comparison to data on neustonic microplastics, due to the inefficiency of benthic and pelagic sampling in collecting such small material in conjunction with the remoteness and size of the benthic and pelagic habitat which may be affected.Furthermore, mechanisms of benthic and belowsurface transport of microplastics are not well understood.Coupling the intrinsic physical properties of microplastics with the ability to simulate hydrodynamic processes in the laboratory and with computer modelling techniques may greatly improve our understanding of the transport pathways of sub-surface microplastics, and likely locations of accumulation.The need for a better means of estimating sub-surface microplastic transport is growing as plastic production levels increase (PlasticsEurope, 2010) and plastic debris accumulation increases worldwide (Barnes, 2009;Wright et al., 2013).Models may help identify and predict regions where ecological communities and fishery-dependent coastal societies are more vulnerable to the potential consequences of plastic pollution, such as associated toxicity to marine organisms and a decline of marine ecological services (Derraik, 2002;Lithner et al., 2009;Lithner et al., 2011).
Approximately half of all manufactured plastics have a higher density than seawater (USEPA, 1992;Morét-Ferguson et al., 2010).HD microplastics are found on beaches, in river and estuary sediments, on continental shelf slopes and in deep-sea benthic environments (Cole et al., 2011).Despite a number of recent studies on benthic plastic pollution, these rarely include information on microplastics due to the difficulties of sampling this size class of material in the deep-sea (Claessens et al., 2011).One extensive study covering European shelf areas reported spatial densities of 0.064-2.63plastic pieces (< 2 cm diameter) per hectare (Galgani et al., 2000).On the California continental shelf, benthic trawls collected microplastics in spatial densities of 6.5 and 1.5 pieces m −3 before and after a storm, respectively (Lattin et al., 2004).Recently, in a study supported by the HERMIONE program, remotely operated vehicle (ROV) video surveys of benthic marine litter in the submarine canyons off the coast of Portugal reported the highest abundances in canyon heads located off the coasts of populated cities (Mordecai et al., 2011).Submarine canyons are known to be conduits for sedimentary material (Monaco et al., 1999;Schmidt et al., 2001;Canals et al., 2006), transporting organic and lithogenic particles from shallow shelf areas to the deep sea via various hydrodynamic processes such as internal waves, tidal circulation, bottom currents, and occasional sediment gravity flows (Van Weering et al., 2002;Canals et al., 2006;de Jesus Mendes and Thomsen, 2007;De Stigter et al., 2007).In light of this, submarine coastal canyon systems may also function as dispersal and transport conduits for benthic microplastics of coastal and fluvial origin (Mordecai et al., 2011).
In this study, we attempt to predict possible microplastic debris transport pathways and likely sites for deposition in a submarine canyon by integrating experimental density, critical shear stress values, and settling velocity properties of preproduction HD plastic pellets into a numerical model.The triple-layer nested model used here was set up with boundary conditions provided by an operational circulation model and an atmospheric forecast model, and utilized a new residence time concept for analysis of the Lagrangian results.Modelling benthic microplastic dispersal from a point source (e.g.river delta, sewage drain, etc.) may be useful in determining the extent and depth to which certain ecosystems are affected, depending on local plastic concentrations, discharge volume, and hydrodynamic conditions.It may also serve political purpose by guiding the development of plastic disposal policy, such as Total Daily Maximum Loads for debris in urban runoff (Moore et al., 2011).

Microplastic hydrodynamic behavior determination
Experimental tests were conducted to determine the average density, settling speed, bed load transport, critical and resuspension shear stress thresholds of HD black preproduction pellets from a sample collected in Los Angeles County, California (received from Algalita Marine Research Institute).
The average spatial dimensions of the pellets were assessed using the software application ImageJ (Rasband, 2012).The average Feret's diameter, defined as the "maximum distance between two points on the selection boundary" (Ferreira and Rasband, 2011, p. 123), of the pellets (N = 350) was measured from photographs modified with a color threshold to allow easy edge determination.The average density of the pellets was determined by measuring the weight of random pellet subsamples (N = 5) and the water displacement of each subsample using distilled water and a graduated cylinder.Average settling velocity of the pellets (N = 50) was determined by video analysis.We extracted a JPG still image each second from the collected video, to allow computation of the velocity of pellets sinking through a 1 m still saltwater column with salinity 36 (PSS), following the method described in Pabortsava et al. (2011).A 20 cm erosion microcosm (Thomsen and Gust, 2000), capable of simulating various benthic shear environments was used to determine the flow velocities at which bed load transport, resuspension, and deposition of pellets (N = 300) occurred.Four replicate runs were conducted according to a predefined calibration table relating rotor angular speed, pump flow and the resultant flow velocity (U * ).Experiments were run in a stepwise manner, in which bottom shear was manually increased over seven, 2-minute duration steps (Table 1).Shear at which 50 % of the particles rolled, slid or saltated was defined as bed load shear velocity U * b .The critical erosion velocity, U * cr , was defined as the shear velocity at which 75 % of the particles were suspended in the water column.U * d , the depositional shear velocity was defined as the flow velocity at which all pellets had settled from suspension.Using water density, the shear velocity values were converted to shear stress values, τ b , τ cr and τ d (N m −2 ) using Eq.(1): where τ (N m −2 ) is shear stress, U * (m s −1 ) is shear velocity, and ρ (kg m −3 ) is the density of seawater (Thomsen et al., 2002).Experiments were filmed to allow for better analysis of particle behavior in laminar flows.

Modelling approach
MOHID Water is a high-resolution numerical model included in the MOHID Water Modelling System, which has been used to model a range of marine environments from coastal areas to open ocean regions (Santos et al., 2002;Braunschweig et al., 2003;Carracedo et al., 2006;Mateus et al., 2012), and which can be used to integrate different oceanographic processes, scales and systems.In this study the model was implemented to simulate the transport patterns of HD black pellets within the Nazaré Canyon.We applied the MOHID hydrodynamic model (Martins et al., 2001;Mateus et al., 2012) to represent the main oceanographic features of the canyon, and coupled it with a Lagrangian model  Braunschweig et al., 2003;Malhadas et al., 2009;Pando et al., 2013) to simulate the transport of the pellets over time.
The model was adapted for the HD black pellets from one used previously in the modelling of organo-mineral aggregate transport (Pando et al., 2013).The hydrodynamic configuration applied in the Nazaré Canyon included three levels of nested models (Fig. 1, Table 2) with one-way coupling as described by Leitão et al. (2005).The three levels each had a vertical resolution of 50 layers; the bottom 43 layers using the Cartesian coordinate system and the top 7 layers using the Sigma coordinate system.The first nested level corresponded to the largest domain, and covered the western Iberian Margin with a spatial grid resolution of 5.6 km.The second level had a spatial grid resolution of 2.0 km and covered the continental shelf area between Figueira da Foz and Ericeira on Portugal's coastline.The third level, and smallest domain, focused specifically on the Nazaré Canyon area at a spatial grid resolution of 400 m.For the model setup, the boundary conditions were provided by the MOHID-PCOMS regional forecasting system (Mateus et al., 2012), forced by the Mercator-Ocean PSY2v2 and linearly composed with the global tide model FES2004 and the atmospheric forecast model MM5 (Dudhia et al., 2005).The physical properties of the operational circulation model were validated on a systematic basis by comparison with satellite remote sensing of the sea surface temperature (SST) data and by deep profiling (Argo) floats data.
To determine whether transport patterns differed by canyon region, the Lagrangian model was used to predict transport of pellets from four 1.44 km 2 monitoring boxes placed within the upper part of the canyon (at depths of 59 m, 262 m, 331 m and 2657 m, box 1-4 respectively) (Fig. 1).Vertical displacement of each box was 0.5 m above the canyon floor and locations of the boxes corresponded to the canyon's head (Box 1), the section of the canyon in line with the shelf break (Box 2), Vitória's tributary (Box 3), and the lower section of the upper canyon (Box 4).Using the results from our experimental investigations on the behavior of the HD black pellets (see Sect. 3.1) and designating the spatial density of the tracers (N = 2000 pellets per box), the model was run for two distinct periods.The model simulated 101 days of the spring/summer period from 1 March 2009 to 24 June 2009 and 101 days of the autumn/winter period from 1 September 2011 to 10 December 2011.The predicted transport and distribution of pellets was computed as the residence time, defined as the temporal interval required by the tracers to leave each monitor box (Pando et al., 2013).

Laboratory experimentation
Average HD black pellet density was 1055 ± 36 kg m −3 (Table 3), approximately 3 % greater than the density of both the seawater (ρ = 1026.69kg m −3 ) used in the laboratory experiments and the seawater (ρ = 1025.1 − 1027.9 kg m −3 ) modelled in the Nazaré Canyon.Accounting for the standard deviation in pellets density, the densest pellets may be up to ∼ 7 % denser than typical seawater (Table 3).Settling velocity of pellets was approximately 28 mm s −1 with little variability between individual pellets (Table 3).HD black pellets displayed uniform resuspension behavior.Bed load transport commenced at shear stresses of 0.014 N m −2 and approximately 50 % of pellets were in bed load transport at ∼ 0.025 N m −2 .75 % of pellets were in suspension at shear stresses of ∼ 0.14 N m −2 and pellets were re-deposited at shear stresses of ∼ 0.087 N m −2 .Shear stress values in Table 3 were approximated from direct observation and video analysis and averaged across replicates.Accuracy of the erosion stress threshold determinations was low due to the slow response time of pellets to changes in flow velocity, and slight differences in individual particle properties (i.e.shape, size, density, degree of bio-fouling).Distinguishing saltation from suspension was also problematic, given the high velocities required to transport the material and the rotational motion of pellets within the chamber.Additionally, the flow of the pump strongly influenced the instantaneous shear within the chamber, with slight fluctuations in pump flow resulting in abrupt changes in transport behavior of the pellets.Maximum shear stress generated in the chamber was limited to ∼ 0.20 N m −2 by the maximum stable pump speed.

Model results
In addition to computed residence time in each monitoring box, three output parameters were used to characterize the pellet transport behavior in the Nazaré Canyon model simulation: distance (the total distance a pellet tracer travelled (km), displacement (the net distance a pellet tracer was transported (km), and velocity (km yr −1 ).The averages of each of  these values, for pellet tracers from each monitor box (during both the spring/summer and autumn/winter modelled periods) are given in Fig. 2. Mean, maximum and standard deviations of modelled values, along with an alternative proxy for overall pellet transport (percentage of pellet tracers which escape each monitor box during modelled period) are given in Table 4.

Residence times of pellets
The residence times of the HD black pellets in each monitor box are depicted in Fig. 3     their monitor box at the culmination of the 101-day model run period.Extrapolating these data, using an average of 95 % of pellet tracers remaining at the end of 101 days allows the average monitor box residence time for microplastics in this study to be estimated at ∼ 5.5 years.The minimal removal of tracers from monitor box 1 indicates that residence times close to shore may be on the decadal scale, assuming hydrodynamic conditions to be similar to those modelled for the entire duration.In contrast, for pellet tracers from monitor box 4 (from which > 5 % of the tracers were modelled to be transported during the 3 month autumn/winter period modelled), the residence time can be estimated to be shorter, at ∼ 3.5 years.Integrating over the entire 200 kmlong canyon system (Tyler et al., 2009), and assuming relatively constant hydrodynamic conditions, the results suggest benthic microplastic transport from the canyon head to the abyssal plain would take place on the centennial or longer scale, if at all.Two signature patterns can be observed in the residence time plots; first, the regular sinusoidal oscillations of the pellet tracers and second, the irregular occurrence of distinct transport events where the fraction of pellets in a monitoring box decreases or increases abruptly (Fig. 3).The regular pellet fraction fluctuations correspond with the diurnal and semi-diurnal component of the tide observed in the Iberian margin region.These sinusoidal fluctuations indicative of tidal forces are observed at all depths for the duration of both simulation periods, but are particularly evident in the shallower monitor boxes 1 and 2 and decrease consistently in amplitude with depth.However, regardless of depth, the tidal water movement appears to have very little effect on modelled net displacement of pellet tracers.Rather, sudden changes (primarily decreases) in the fraction of tracers in monitor boxes suggest transport occurs primarily within the canyon during occasional periods of elevated oceanographic flow conditions, e.g.sediment gravity flows or internal waves which are reported to occur in the canyon (De Stigter et al., 2007;Quaresma et al., 2007;Muacho et al., 2013).These sudden transport events often occurred simultaneously in each box, signifying that these pellet tracer movements were likely the result of a large-scale event, rather than small-scale surface disturbances.
One exceptionally strong transport event can be seen in Fig. 3b, where 10 % of pellets were transported out of box 2 on 1 October 2011 of the autumn/winter model run.To analyze the cause of this strong dispersion event, which was identified to occur between 18:00 and 19:00 UTC, we www.biogeosciences.net/10/7957/2013/Biogeosciences, 10, 7957-7970, 2013  investigated the model input data for several days preceding the event.As wind shear stress could influence pellet transport at 262 m depth, wind patterns were analyzed but these did not show any significant change in direction or intensity.River discharges have been reported to influence current regimes on the western Iberian margin (Oliveira et al., 2007;Martín et al., 2011); however, flow data provided by the hydrometric stations (www.snirh.pt)from the Douro, Mondego and Tagus rivers near the canyon head showed that the contribution from those rivers' discharges was not significant during this period.Lastly, the modelled bottom currents in box 2 were analyzed and these showed a consistent pattern in the area during this period, with peak horizontal and vertical velocities in the bottom layer.The horizontal and vertical velocity modulus along the canyon axis for the time period between 17:00 and 20:00 UTC on 1 October 2011 is plotted in Fig. 4. Velocities up to 0.3 m s −1 were modelled to occur in the bottom layer in box 2, which are sufficient to suspend and transport pellets.Bottom shear velocities near box 2 during this time period were a factor of 3 times higher than those used to suspend HD black pellets during laboratory experimentation (Tables 1 and 3).Considering water density (ρ = 1026.69kg m −3 ) used in laboratory simulations are similar to modelled water density at box 2 (σ T ∼ 27 kg m −3 ) the shear stresses experienced by the pellets in the model should be comparable to laboratory experimentation.To determine whether or not the increased current velocity was an artifact of an internal wave, we plotted the isopycnal surfaces, σ T (kg m −3 ), throughout the water column and along the axis of the canyon, arraying a series of plots for the same time period on 1 October 2011 (Fig. 5).The propagation of an internal wave with amplitude up to 200 m is visible approaching and breaking within the upper canyon.To further verify that the pellet transport was caused by the passing of the internal wave, we compared the mean isopycnal surfaces and mean velocity modulus along the canyon axis between the 1 October and a day during which pellet transport was modelled not to occur: 19 October 2011 (Fig. 6).The increased bottom shear stresses near box 2 at the shelf break coincide with the development of an internal wave at the head of the canyon (Fig. s 4 and 5).Investigation of the bottom current velocity and isopycnal surfaces on the 1 and 19 October 2011 suggest that mean current velocities are insufficient to induce unidirectional transport of pellets, but amplified current velocities due to large-scale hydrodynamic events, such as internal waves, are sufficient.

Seasonal variability
As mentioned in Sect.3.1, abrupt displacements of pellet tracers from monitoring boxes were more evident in the residence time plots (Fig. 3) during the autumn/winter model run than in the spring/summer run.The other modelled transport parameters (distance, displacement and velocity) also differed by season, particularly for the shallow monitor boxes.
Distance and displacement values were consistently higher during autumn/winter, up to factors of 5 or 6 times higher for pellet tracers originating in boxes 1 and 2 (Fig. 2).The standard deviations in these values were consistently higher during the autumn/winter period, as were the maximum modelled values, a trend indicative of a more chaotic and variable hydrodynamic environment (Table 4).Average pellet tracer velocities ranged between ∼ 0.1 and 1.0 km yr −1 in the spring/summer period and ∼ 0.2 to 0.9 km yr −1 in the autumn/winter period suggesting that pellet velocity did not change significantly between seasons (Table 4).However, the average tracer velocity in the autumn/winter (0.58 km yr −1 ) was slightly greater than in the spring/summer (0.44 km yr −1 ).As with the other modelled parameters, the maximum values and standard deviations making up this average were also higher during the autumn/winter period (Table 4).
In the model runs, large-scale transport events appear to occur more regularly and with similar impact in each monitor box in the spring/summer period as compared to the autumn/winter period.During the autumn/winter, when these sudden events were modelled to occur, they impacted pellet transport to different degrees in each monitoring box.For example, the large transport event modelled to occur on day 30 of the autumn/winter run effectively removed ∼ 10 % of pellet tracers from box 2, but had only a minimal effect on concentrations within boxes 1, 3 and 4 (Fig. 3).

Differences in transport by canyon location
From the model output during both spring/summer and autumn/winter periods, pellet tracers throughout the canyon travelled greater distances than they were displaced (Fig. 2), further indicating the tracers were transported in an oscillating manner, up and down canyon repeatedly, as suggested by the residence times shown in Fig. 3.
During the spring/summer period, boxes 1 and 3 exhibited similarly low displacements, distances and velocities of pellet tracer transport (Fig. 2), but had the largest difference in the percentage of tracers that ultimately escaped the monitor box (Table 4).Similarly, transport behavior of pellet tracers in box 2 and box 4 was similar, but the percentage of tracers escaping from monitor box 4 was a factor of 3 times higher than that modelled for box 2. These incongruences suggest that pellet movement within and pellet export from the monitor boxes are not necessarily correlated.Alternatively to the spring/summer period, during the autumn/winter period, pellet tracers in boxes 1, 2, and 4 exhibited consistently high distance, displacement and velocity behaviors as compared to pellet tracers in box 3 (Fig. 2) but overall, the percentage of escaped tracers increased with depth in the autumn/winter modelled run, and, with the exception of box 3, in the spring/summer modelled run (Fig. 3, Table 4).

Canyon topography and flow
Cross sections through the canyon within each monitoring box, with average modelled flow velocities for 28 May 2009, are given in Fig. 7.As the figure shows, much of the modelled high velocity flow takes place in the shallow waters in monitoring boxes 1 and 2 for the day presented.Where the canyon starts to open up, near monitor box 3, higher velocity flows are modelled to occur at greater depths too.In box 4, where the canyon opens up considerably, chaotic higher velocity flows can be observed at various depths, particularly at around 1000 m.In boxes 1 and 2, the cross canyon flow evident in the surface waters does not influence bottom flow, possibly constrained by the deep, narrow canyon topography.Alternatively, in boxes 3 and 4 where the canyon topography is of lower relief, near-surface flow patterns can propagate to deeper waters (in box 3, Fig. 7c particularly), potentially increasing the frequency of bottom velocities sufficiently high for pellet resuspension.

Discussion
The differences in pellet transport (both within and between seasons) modelled for each monitor box location reflect the changes in the hydrographic regime over time and the topography within the Nazaré Canyon.
The model simulations presented here indicate that the dispersion of ∼ 5 mm diameter HD microplastic pellets through the Nazaré Canyon during typical spring/summer season conditions is likely slow, but may increase with the Fig. 8.The critical bed shear stress erosion curve for quartz relates particle (sediment) size to critical shear stress, τ cr , and includes average diameter (d 50 ) 4000 µm benthic boundary layer aggregate data point (Thomsen et al., 2002).The mean HD black pellet size (d 50 ∼ 4700 µm) and τ cr is plotted over the curve for comparison of aggregate and plastic erosional behavior.
intensification of the hydrographic regime during the autumn/winter season.The model output suggests that topography restricts pellet movement at the head of the canyon making it a potential accumulation area for non-buoyant plastic debris, but allows escaped pellets to disperse more quickly at the shelf break and in deeper sections with the occurrence of large-scale hydrodynamic events due to the widening and deepening of the canyon axis.The canyon, situated on the western Iberian Margin of the eastern North Atlantic Ocean, is a hydrographic region characterized mainly by tidal currents, internal waves, and upwelling (Vitorino et al., 2002;Quaresma et al., 2007).Throughout the upper canyon, pellet movement appears from the model to be consistently affected by tidal forces as can be seen in Fig. 3; with predicted residence time of pellets fluctuating in a sinusoidal pattern characteristic of the peaks of the M 2 tide.Recirculation of tidal currents within the canyon has been shown to actively resuspend organic and fine-grained lithogenic material (de Jesus Mendes and Thomsen, 2007;De Stigter et al., 2007).From our model, it appears that tidal bottom shear stresses may also be sufficiently high, not for resuspension, but for bed load transport of HD plastic pellets in the upper canyon and to a lesser degree in the lower canyon.Near the shelf break, (boxes 2 and 3 at ∼ 300 m) where the canyon axis begins to deepen and widen (Fig. 7), the bottleneck shape (Fig. 1) of the canyon may be a source of higher current velocities and wave-induced turbulence, resulting in increased pellet transport, down-canyon.Tidally induced internal waves moving shoreward may amplify as they move into the shal-low, narrow canyon head, causing them to destabilize at the shelf break (Fig. 4), generating turbulence and elevated bottom shear stresses sufficient to resuspend microplastic debris on the seafloor.The regular and synchronized occurrence of sudden transport events in the three deepest monitoring boxes suggests that these events correlate with a regular hydrodynamic event, such as tidally induced internal waves, with associated sediment gravity flows and resuspension of bottom material.The variable effect of such events on pellet transport in each monitor box indicates that the local canyon topography depth also affects the degree to which the largescale hydrodynamic disturbance can transport pellets downcanyon.Other similar phenomena, such as dense water cascades and turbidity flows, which can occur within canyon systems due to large-scale hydrodynamic processes, are not considered here as they are not known to occur in the Nazaré Canyon on yearly timescales (De Stigter et al., 2007).The model run results support the hypothesis that in certain areas of the canyon internal waves are likely significant contributors to resuspension and mobilization of deposited materials, such as microplastics.
Sedimentology within the Nazaré Canyon, and specifically transport of organo-mineral aggregates (OMAs) of various size classes, has been thoroughly investigated (Schmidt et al., 2001;Thomsen et al., 2002;Van Weering et al., 2002;De Stigter et al., 2007;Oliveira et al., 2007;Arzola et al., 2008;Lastras et al., 2009;Pando et al., 2013).Comparing simulated pellet transport to OMA transport in a comparable model study in the Nazaré Canyon (Pando et al., 2013), indicates that microplastic pellets were modelled to behave more similarly to OMAs of size classes 2 mm and 4 mm than those of 429 µm.Modelled transport predictions also indicated that pellets were transported for shorter distances than OMA aggregates, likely due to the considerably higher settling velocity and critical erosion shear stress values of the pellets.In Table 3, the physical parameters of comparably sized OMAs from the Iberian continental margin are listed for comparison to HD black pellets (Thomsen et al., 2002;de Jesus Mendes and Thomsen, 2007).The pellets have a settling velocity 7 times greater and erosional shear stresses approximately 5 times greater than OMAs.In Fig. 8, the erosional shear stress of the pellets and OMAs are plotted on a quartz erosion curve as taken from Thomsen et al. (2002).The position of the pellets on the curve indicates that they can be expected to behave more similarly to compact, high-density medium-grain sand particles (d 50 ∼ 600 µm) than to loosely packed, low-density benthic boundary layer aggregates of 4 mm size when under the influence of bottom currents.This is in agreement with the model results, where average pellet velocities are < 1.0 km per year.Given that the length of the canyon is ∼ 200 km (Tyler et al., 2009), transport of benthic microplastics from the shelf area to the abyssal plain would not occur during the timescales modelled here, unless greatly enhanced by large sediment gravity flow events.

A. Ballent et al.: Modelled transport of benthic marine microplastic pollution
During autumn and winter seasons, storms may cause sporadic sediment gravity flow events and higher river discharges (Martín et al., 2011), neither of which were specifically included in the model parameters.Such events, observed to occur on yearly scales, may transport high volumes of fine-grained sediment and debris from the upper to the middle and lower canyon.Sediment gravity flows strong enough to transport sandy sediment likely only occur only on centennial timescales (De Stigter et al., 2007), and following such events, plastics not wholly transported from the canyon would likely be buried within deposited sediment, further increasing residence time in the canyon system (Galgani et al., 1996;Mordecai et al., 2011).The rarity of flows of such magnitudes suggests that microplastics will accumulate within the upper canyon, assuming that the majority of plastic waste comes from land, a hypothesis supported by in situ data collected by Mordecai et al. (2011), wherein a correlation between macro-debris and distance from the coastline in the Lisbon and Setubal submarine canyons is reported.
Benthic debris may also accumulate in certain zones of the canyon where further transport is inhibited by benthic topography, as suggested by Galgani et al. (1996) and Mordecai et al. (2011).Models could be used to locate and identify such depositional areas given that high-resolution physical oceanography and topography data is available.Additional field data (e.g.plastic counts from box core sediment samples and near bottom sediment trap samples) should be used to validate dispersal and depositional predictions.
Critical erosion values in this model were determined in a laminar flow environment by simulating the benthic boundary layer velocities found in deep sea environments (Thomsen et al., 2002).However, flows generated by tides, waves and uneven bottom surfaces may result in turbulent conditions (De Stigter et al., 2007;Martín et al., 2011), and consequently a thinner benthic boundary layer and changes in the resuspension behavior of microplastics.Wave action, tides, turbidity flows and storm events may all induce turbulent bottom shears, particularly in shallow areas of a canyon system.Future research should incorporate both laminarand turbulence-induced critical erosion shears, in connection with location and topography, to improve the accuracy of future model predictions of benthic plastic transport.

Conclusions
This investigation was an attempt to gauge the degree to which the intrinsic properties of plastic debris affect their transport within the sub-surface waters of a modelled marine canyon environment.It is crucial to understand how microplastics are transported below the surface to more accurately estimate global distribution, residence times, convergence zones and ecological consequences, which can be done using numerical models.The model presented here indicates slow transport of benthic microplastics in the Nazaré Canyon, implying long-term exposure of benthic canyon ecosystems to plastics, but also suggests that large-scale water movement, such as those associated with storms and internal waves, may intensify down-canyon transport in episodic short duration events.Due to the long residence times for benthic microplastics indicated by our study, future investigations into the rates at which benthic microplastics degrade are important to better quantify realistic residence times.Potential changes in the intrinsic properties of microplastics as they degrade toward nanometer scale, may alter transport properties and, in turn, the results of long run-time transport models.Future research should also focus on the ecological consequences of plastic exposure in benthic environments, particularly in critical areas such as biodiversity hotspots, to allow the development of preventative measures and policy/legislation changes if required.Decreasing the amount of plastic debris originating from urban consumers would reduce exposure levels in many deep sea regions close to shore, such as the canyon ecosystems focused on in the current study.
With the culmination of this study, it is clear that further research is needed to accurately estimate the amount of plastic residing in the benthic oceans and to understand the sub-surface behavior of non-buoyant microplastics and their sinks within the natural environment.

Fig. 1 .
Fig. 1.The nested domains of the hydrodynamic model used in the simulation of the dispersal of plastic preproduction pellets in the Nazaré Canyon on the western Iberian Margin: (a) first level: MOHID-PCOMS; (b) second level: Figueira da Foz-Peniche; (c) third level: Nazaré Canyon and the locations of the monitor boxes 1-4 at depths 59 m, 262 m, 331 m and 2657 m along the canyon axis.
as the fraction of pellet tracers remaining inside the monitor box over time.In all cases of both the spring/summer and autumn/winter model runs, ∼ 90-100 % of the pellet tracers were not transported outside Biogeosciences, 10, 7957-7970, 2013 www.biogeosciences.net/10/7957/2013/

Fig. 3 .
Fig. 3.The residence times of the HD black pellet tracers in each monitor box of the Nazaré Canyon model simulation in (a) the spring/summer 2009 model run and (b) the autumn/winter 2011 model run.

Fig. 4 .
Fig. 4. The velocity modulus * vertical currents (m s −1 ) are plotted hourly from 17:00 to 20:00 UTC (a-d, respectively) for 1 October 2011, along the Nazaré Canyon axis (red curve on top) during which a large export of pellet tracers was observed in monitor box 2.

Fig. 5 .
Fig. 5.The propagation of a modelled internal wave hourly from 17:00 to 20:00 UTC (a-d, respectively) on 1 October 2011 in the Nazaré Canyon where the vertical distribution of water density, σ T (kg m −3 ), at depths up to 200 m varies rapidly along the canyon axis, particularly at the canyon head.The contour interval for the isopycnal surfaces (black lines) is 0.3 kg m −3 except for the range (27.0-27.1)where the step is 0.5 kg m −3 .

Fig. 6 .
Fig. 6.Comparison of vertical distribution of water density σ T (kg m −3 ) and velocity modulus * vertical current (m s −1 ) averages for two days, 19 October (a and c) and 1 October (b and d) 2011, along the Nazaré Canyon axis (red curve on top).In (a) and (b), the contour interval for the isopycnal surfaces (black lines) is 0.3 kg m −3 except for the range (27.0-27.1)where the step is 0.5 kg m −3 .

Fig. 7 .
Fig. 7. Average modelled flow data, velocity modulus * vertical currents (m s −1 ), for 28 May 2009, is plotted for the cross section through the canyon system at each monitor box location.

Table 1 .
Resuspension steps of shear velocity for bed load transport, resuspension and deposition within erosion microcosm containing saltwater of density ρ = 1026.20kg m −3 .

Table 2 .
Main characteristics of the hydrodynamic nested models including coordinates, dimension, resolution and time step used in simulation of HD black pellets in the Nazaré Canyon.

Table 3 .
(Thomsen et al., 2002;de Jesus Mendes and Thomsen, 2007)ts from the Los Angeles beach sample.Settling velocity, bed loadτ b , criticalτ cr and depositionalτ d shear stress of the HD black pellets as tested in saltwater of ρ = 1026.69kgm−3 .Mean diameterd 50 , density, settling velocity, critical and depositional shear stresses of BBL aggregates(Thomsen et al., 2002;de Jesus Mendes and Thomsen, 2007)for comparison to HD black pellets.Standard deviations indicate differences between measurement replicates.

Table 4 .
Average distance, average displacement and average velocity of pellet tracers and percentage of pellet tracers transported from the monitor boxes by the end of each Nazaré Canyon model simulation (S/S: Spring/Summer, A/W: Autumn/Winter).N = 2000 pellet tracers per monitor box.