Modelling size distributions of marine plastics under the influence of continuous cascading fragmentation

Field studies have shown that plastic fragments make up the majority of plastic pollution in the oceans in terms of abundance. How quickly environmental plastics fragment is not well understood, however. Here, we study this fragmentation process by considering a model which captures continuous fragmentation of particles over time in a cascading fashion. With this cascading fragmentation model we simulate particle size distributions (PSDs), specifying the abundance or mass of particles for different size classes. The fragmentation model is coupled to an environmental box model, simulating the distributions of plastic particles in the ocean, coastal waters, and on the beach. We compare the modelled PSDs to available observations, and use the results to illustrate the effect of size-selective processes such as vertical mixing in the water column and resuspension of particles from the beach into coastal waters. The model quantifies the role of fragmentation on the marine plastic mass budget: while fragmentation is a major source of (secondary) plastic particles in terms of abundance, it seems to have a minor effect on the total mass of particles larger than 0.1 mm. Future comparison to observed PSD data allow us to understand size-selective plastic transport in the environment, and potentially inform us on plastic longevity.


Introduction
Studies have shown that fragments make up the majority of marine plastic litter in terms of abundance (Cózar et al., 2015;Suaria et al., 2016). The large amount of fragments is evident from particle size distribution (PSD) data, specifying the abundance or mass of particles for different size classes. An overview of PSD data from various studies is given in Kooi and Koelmans (2019); some examples are presented in figure 1. What is commonly observed in PSD data is a power law for larger fragments (>1 mm in figure 1), resulting in a straight line on a log-log scale. Oftentimes, a maximum in the PSD is observed at smaller particle sizes (∼1 mm in figure 1), ending the power law regime.
If we want to explain the particular shapes of measured PSD data, we should first investigate the fragmentation process. Fragmentation of plastics is likely dominant on beaches, where plastics are subjected to UV-radiation, oxidation, and higher temperatures, embrittling the particles, which enhances the breaking down of particles by mechanical abrasion (Andrady, 2011;Kalogerakis et al., 2017;Song et al., 2017;Efimova et al., 2018). Fragmentation models have been proposed in e.g. Cózar et al. (2014), hypothesising that the PSD slope depends on whether particles break down in a three-dimensional fashion (i.e. like a cube), or more in a two-dimensional fashion (like a thin sheet).
While the main driver behind the PSD might be fragmentation, physical processes can have a size-selective influence on plastic particles . Vertical mixing has been shown to influence PSDs measured in the water (Reisser et al., 2015;Poulain et al., 2019). Assuming similar particle shapes and densities, smaller particles tend to have lower rise velocities due to a lower buoyancy which counteracts the water friction. Turbulent mixing, induced by for example the wind, will therefore distribute these smaller particles across larger depths (Kukulka et al., 2012;Chor et al., 2018), reducing smaller size fractions in PSDs measured by nets at the ocean surface (typically submerged ±10-50 centimeters depending on net type, see e.g. Pedrotti et al. (2016), Cózar et al. (2015), and Suaria et al. (2016)).
Size-selective lateral transport of plastic particles is likely to influence the PSD shape as well. Floating particles experience a net drift caused by waves at the ocean surface, i.e. Stokes drift (Stokes, 1847). Bigger (more buoyant) particles likely experience more influence from Stokes drift, given its limited depth of influence (Bremer and Breivik, 2017;Breivik et al., 2016). Model studies have shown that Stokes drift tends to push particles towards coastal areas (Iwasaki et al., 2017;Onink et al., 2019;. In Isobe et al. (2014) it was shown that PSD data obtained close to the coast indeed showed relative overabundances of larger plastic particles. In figure 1 this is illustrated by differentiating between coastal and offshore areas: relatively high abundances of small fragments in the open ocean can be seen.
Coastal processes, such as beaching and resuspension, can be size-selective. In Hinata et al. (2017), residence times of particles on beaches were estimated using tagged litter. Higher particle rise velocities in the water were related to longer residence times, as these particles are more likely to be pushed to the backshore by wave swash. This could mean that larger objects remain longer on beaches, and hence experience more weathering (Hinata et al., 2020).
Finally, PSDs could be influenced by size selective sinking, induced by for example biofouling (Ryan, 2015). Biofouling models predict that smaller particles, which have a larger surface to volume ratio, tend to sink more quickly . This has been observed in experimental studies as well (Fazey and Ryan, 2016).
In this work, we first consider a fragmentation model, based on fractal theory (Turcotte, 1986;Charalambous, 2015). We will use this model to simulate how PSDs resulting from fragmentation might evolve over time. Afterwards, we couple this fragmentation model to a conceptual box model, simulating movement   Cózar et al. (2014Cózar et al. ( , 2015, Ruiz-Orejón et al. (2018), and Isobe et al. (2014Isobe et al. ( , 2015. Coastal samples are defined as less than 15 kilometers from the shoreline. All particle size distributions have been normalized relative to the total abundance of large items (> 10 mm) to show differences for small particle sizes clearly. Most datasets seem to follow a power law for the larger particle size classes (see figure 7 below). Coastal samples tend to have relatively few small particles (or: relatively many large particles), and all distributions show a peak in abundance around 1 mm, instead of a monotonic increase with smaller sizes. of plastic particles between the beach, coastal water, and open ocean. Different particle sizes experience different forcings in this model, under the influence of wind-induced turbulence (Kukulka et al., 2012;Poulain et al., 2019), sizeselective transport due to Stokes drift (Breivik et al., 2016), and size-selective resuspension from beaches (Hinata et al., 2017). We will compare the modelled PSDs to available observations. Instead of counting the amount of particles in a given size-range, one can also measure their mass. In this work, we make a distinction between the two, by considering a number (i.e. abundance) size distribution (NSD), and a mass size distribution (MSD). The term PSD is used when applicable to both.

The cascading fragmentation model
The fragmentation model discussed here is based on simple fractal geometries. We define the spatial dimension as D N . When D N = 3, we start with a cube with a size of L × L × L which we call the parent object, see figure 2a. This cube can be split in eight equally-sized cubes, which can each be recursively split again. The size class of the parent object is defined as k = 0, the size class of the cubes with length L/2 is defined as k = 1, and so on. When D N = 2, the starting object is a sheet instead of a cube, which can be split in four smaller Figure 2: a) Illustration of a parent cube (size class k = 0), consisting of successively smaller cubes, based on Turcotte (1986). Only 3 iterations until size class k = 3 are shown here, the size class can increase indefinitely. b) Illustration of the cascading fragmentation model with p = 0.5 and after one fragmentation event i f , based on Turcotte (1986) and Charalambous (2015) sheets each time the size class increases. Cózar et al. (2014) presented a fragmentation model where objects are broken down into a set of smaller (equally sized) fragments in a series of successive fragmentation events. This fragmentation model was used to explain why measured PSDs often resemble power laws, i.e. functions of the form where n is the abundance, l is the particle size, α is the power law slope, and C is a constant. However, this fragmentation model requires a constant input of new parent objects to achieve a power law, while laboratory experiments have shown that power laws in the PSD also appear after fragmenting a single input of parent objects (Song et al., 2017). The fragmentation model used here builds upon the work of Turcotte (1986), where it was noted that scale-invariance of the fragmentation process, whether it be caused by weathering, explosions, or impacts, leads to such a power law. The idea behind the model of Turcotte (1986) can be illustrated using figure 2b. Following one fragmentation event i f , a certain fraction p of the original cube (size class k = 0) splits off. For example, if p = 0.5, this results in 4 fragments of k = 1 splitting off, leaving 0.5 object in size class k = 0. This process is assumed to be scale invariant: a fraction p will split off from the fragments in size class k = 1 as well: 16 fragments of size class k = 2 are created, and 4 × 0.5 = 2 fragments are left in size class k = 1. This process is repeated indefinitely. Bird et al. (2009) and Gregory et al. (2012) extended this model by including a temporal component, with each fragment breaking down further as i f progresses. Charalambous (2015) showed that repeatedly breaking down fragments over discrete steps of i f is a sequence of independent and identical Bernoulli trials with a chance of success p, yielding a negative binomial distribution. This is rewritten in terms of a continuous fragmentation index f (instead of the discrete i f ), yielding a probability density function giving the mass m in size class k at fragmentation index f as: where Γ is the gamma function. We will call this the cascading fragmentation model. We assume that f is directly proportional to time in the environment, and will review this assumption in the discussion. We assume p to be a material property. The amount of fragments in a given size class is estimated by multiplying the mass with 2 D N k , a factor determining how many fragments of size class k fit inside the parent object: We use D N = 3 as the baseline. However, this factor is D N = 2 for purely flat objects like plastic sheets and D N = 1 for fibers or lines. As real-world samples contain a combination of these objects, the value for D N in the environment can be a non-integer between 1 and 3. Figure 3a shows the NSD resulting from the cascading fragmentation model at various fragmentation indices f . We start with one cube with a length L of 1 mm at f = 0. The continuous description in (3) allows us to model the amount of fragments at a very small fragmentation index of f = 0.001. There are few larger fragments (> 10 −2 mm) per parent object at this stage. At f = 1 we have exactly a power law in the NSD, equivalent to the model by Turcotte (1986). A fractal dimension D f of the object formed by all fractions can be defined, relating to D N and p by: The NSD power law slope at f = 1 is given by this fractal dimension. Fragments can be broken down further, eventually resulting in the NSD shown for f = 10. This is not a power law anymore, and the slope of this curve has increased significantly, with relatively many particles in the small size classes. The NSD (units: n) can be normalized, by dividing the amount of fragments by the size class bin width (units: n mm −1 ). These normalized NSDs are presented by the dashed lines. Because of the log-scale on the x-axis, the distance between the given particle sizes increases by a constant factor. This increases the magnitude of the normalized NSD slopes by 1 compared to the discrete NSD. These normalized NSDs allow for comparison between different studies, and are therefore usually reported in literature. Figure 3b shows the same analysis in terms of mass, i.e. the MSD, starting with one cube of 1 g and 1 mm 3 . As fragmentation progresses, mass shifts from the large fragments towards smaller fragments. At f = 1 we have a power law: the difference in the slope between the NSD and MSD is 3, resulting from the 2 D N k term in (3), with D N = 3.

Environmental box model
With the cascading fragmentation model we can now simulate PSDs over time. How particles are affected by the environment varies for different particle sizes, however. The combination of fragmentation and size-selective transport is investigated using a box model, presented in figure 4. The boxes in this model represent three different environmental regions: the beach, coastal waters, and open ocean.
Particles can move between the different environments, defined by a set of transition probabilities (P ): particles can move to a different environmental box (the arrows on the right in figure 4), remain in the current box (recurring arrows on the left), or vanish from the system (i.e. a sink, arrows on the left). Subscripts in figure 4 denote ocean, coast, beach, or sink (O, C, B, and S respectively).
Besides different environmental regions, we have different particle size classes. For a given size class, certain mass fractions will move to smaller size classes under the influence of fragmentation. These fractions are estimated by evaluating (2) for the given time step of the box model. Similarly, (3) is evaluated to determine the abundance of fragments moving to smaller size classes. Fragmentation is assumed to only happen on the beach, where degradation is expected to be much more effective than in the sea (Efimova et al., 2018). Particles move between and within the ocean, coast, and beach boxes, and for each box there is a probability of being removed from the system.
Each environmental box contains a range of different particle size classes. The combination of environmental transition probabilities and fragmentation is modelled using a transition matrix. For example, taking 15 different size classes and 3 environmental compartments leads to a transition matrix of size 45 × 45. Further details are given in the supplementary material (section S2).
Environmental transition probabilities are based on previous findings from literature, and are discussed in the next paragraphs. A set of baseline parameters are varied in the box model to investigate their influence on the modelled PSDs. Where possible, box model parameters are calibrated for the Mediterranean sea, for which plastic sources and transport estimates are available (Kaandorp et al., 2020). For fragmentation we have two unknown parameters: p in (2), and a fragmentation rate λ relating f to time in the environment. These are estimated by comparing the model to experimental weathering data from Song et al. (2017), presented in the results section.
We assume that new plastic objects are introduced on the beach, and assume an initial length of 200 mm based on typical dimensions of municipal plastic waste, see the supplementary material (section S4). We use D N = 3 as the baseline, i.e. cubical-shaped objects. The model time step is set to one week.
Ocean Transport A Lagrangian simulation of floating plastic in the Mediterranean, see Kaandorp et al. (2020), is used to determine transition probabilities within and between the ocean and coast (P O,O , P O,C , P C,C , and P C,O ). The coast is defined as the ocean within 15 kilometers of the coastline.
The baseline transition probabilities are set to a uniform value for all parti-cle sizes, assuming they reside at the ocean surface where they experience the maximum Stokes drift. We will then evaluate the influence of size-selective lateral transport, induced by vertical mixing of particles and differences in Stokes drift over depth, see Isobe et al. (2014) and Iwasaki et al. (2017) for similar approaches. Vertical mixing of particles is estimated using the approach from Poulain et al. (2019): for given particle sizes, the rise velocity w b is estimated by assuming ellipsoid particles with a density of 950 kg m −3 . From this rise velocity we estimate the median particle depth, using the particle density profiles from Kukulka et al. (2012). The Stokes drift is estimated at this depth, assuming a Stokes profile based on the Phillips wave spectrum, see Breivik et al. (2016). Beaching and resuspension Using the same approach as in Kaandorp et al.
(2020), we estimate P C,B by analysing drifter buoy data: from a set of 1682 drifters in the Mediterranean (Menna et al., 2017), we calculate how much time these drifters spend near the coast before beaching. For drifter buoys within 15 kilometers of the coastline, the beaching rate is estimated to be about 6.69 · 10 −3 day −1 (corresponding to an e-folding time-scale τ CB of 149 days). In Kaandorp et al. (2020) it was estimated that τ CB for plastic particles is about 3 times lower than that for drifter buoys. We will therefore use τ CB = 50 days as the baseline estimate here. In Hinata et al. (2017), tagged litter on beaches was tracked over 1-2 years, to estimate their residence times τ BC . As a baseline estimate we use τ BC = 211 days, reported for small plastic floats (corresponding to P BC = 4.73· 10 −3 day −1 ). We will also look at the effect of size-selective resuspension, for which the empirical relation from Hinata et al. (2017) is used, i.e., τ BC (w b ) = 2.6 · 10 2 w b + 7.1, where τ BC is given in days, and w b in m s −1 .
Sinks The box model also requires transition probabilities for removal of particles: P O,S , P C,S , P B,S . We assume these are the same in all compartments, denoted by P S . A given value for P S yields a certain amount of steady-state mass in the system. We take the estimated input of waste into the Mediterranean from Kaandorp et al. (2020) (2,500 metric tonnes for the year 2015), and the estimated total floating mass from Cózar et al. (2015) (2,000 metric tonnes). The value for P S is iterated until this mass balance is satisfied, see the supplementary material (section S2) for more information.

Calibration with laboratory experiments
In Song et al. (2017), plastic pellets were subjected to different levels of UV exposure and to 2 months of mechanical abrasion with sand, simulating a beach environment. We use these data to calibrate the cascading fragmentation model. The data for polyethylene (PE) and polypropylene (PP) pellets (26 mm 3 and 19 mm 3 respectively) are used, as these are the most abundant polymers in surface ocean (Bond et al., 2018).
We assume a single p value per material, and D N = 3. The fragmentation index f is allowed to vary between the different levels of UV exposure when fitting the data. By fixing p and varying f , we get a robust estimate for the unknown parameter p for which we need a plausible value in the box model. We can expect that f is larger for particles subjected to longer periods of UV exposure, since embrittlement will make it easier for the mechanical abrasion to wear down the particles.
Resulting NSD fits using weighted least squares are presented in figure 5, top row, fitted values for f are presented in the legend. For PE particles, the best fit results in p = 0.39, for PP particles p = 0.45. The experimental data are still at the early stage of fragmentation (f < 1), with few fragments per parent pellet, except for small fragment sizes.
There is a good fit for the PE data, with almost all simulations within the data errorbars (one relative standard deviation). For PP there is a good fit for 0, 2 and 6 months of UV exposure. At 12 months of UV exposure there is more mismatch for the smallest size class (0.05-0.10 mm). This is also the only case where the estimated f is lower than for the previous level of UV exposure.
The bottom row of figure 5 compares the estimated volume fractions of the parent pellets and the fragments. Generally, the modelled volume fraction of the parent pellet is estimated well, although there is some overprediction for PE with 12 months of UV exposure. The modelled fragment volumes are higher than the ones estimated in Song et al. (2017). A possible explanation is that some of the larger fragments, which contribute to most of the volume when f < 1 (see figure 3), could have been missed in the experimental setting since there are very few of these per parent object (e.g. tenths or hundredths).
Following this analysis, we set p to 0.4 in the box model, within the range of fitted p for PE and PP. A fragmentation rate λ needs to be chosen, specifying how much f increases per unit time. Following Song et al. (2017) we assume that 12 months of laboratory UV exposure is roughly related to 4.2 years of environmental exposure. Taking our estimated fragmentation indices for PE and PP results in fragmentation rates λ of 1.8 · 10 −2 f year −1 to 6.9 · 10 −2 f year −1 . The value for PE is used as the baseline here, given it is the most common environmental polymer (Bond et al., 2018). We acknowledge that this fragmentation rate is still very uncertain, and f might vary non-linearly with time in reality.

Modelled environmental particle size distributions
Now that we have estimates for the transition probabilities in the box model and estimates for the fragmentation parameters, we can simulate PSDs in the environment. We will quantify the power law slope α of the results by numerically maximizing the log-likelihood of the data (Virkar and Clauset, 2014): where b are the bin boundaries used to discretize the data, containing n i samples in the bin with index i, and n = n i . In some cases, not the entire particle size range adheres to a power law. The lower bound of the power law domain is estimated by minimizing the Kolmogorov-Smirnov statistic between the modelled NSD and the theoretical power law NSD (Virkar and Clauset, 2014).
NSDs resulting from the box model are shown in figure 6, corresponding MSDs and a table with parameter settings can be found in the supplementary material (section S1). Fragmentation is expected to increase the fraction of small particles, increasing α over time (see figure 3). However, environmental sinks limit the magnitude of α: assuming a constant removal rate of plastic particles, smaller fragments, which tend to be older, have a higher probability of being removed from the environment. This combination of fragmentation and environmental sinks eventually leads to an equilibrium, or statistical steady state. This is illustrated in figure 6a using the box model with the baseline parameters described in section 2.2. As time progresses, the relative proportion of fragments to parent objects increases. In this scenario, it takes on the order of years for the NSD to resemble the steady state (red dashed line). The magnitude of the environmental sinks is high enough to avoid long persistence of fragmented particles: there are still relatively many parent objects, and α = 2.57 is still below the value derived from the fractal dimension of α = 2.67 from (4).
Steady state NSDs for different scenarios are presented in figure 6b. Results for the baseline parameters (blue lines) almost overlap with the results where size-selective lateral transport is added to the box model, induced by vertical mixing and Stokes drift (orange lines). In the baseline scenario α = 2.57 for all three NSDs. When adding size-selective ocean transport, larger particles tend to move more frequently from the ocean to the coast. This results in slightly more small particles in the ocean box, increasing the power law slope here to α = 2.73. Adding size-selective resuspension of particles (Hinata et al., 2017) has a strong effect (green lines). Bigger objects have longer residence times on the beach, and therefore undergo more fragmentation. This produces a large number of smaller fragments with shorter residence times, which therefore move more rapidly to the coastal and ocean cells. This near-shore trapping of larger plastic objects was already hypothesized in e.g. Isobe et al. (2014). The empirical resuspension relation (5) causes the model to deviate from a power law, the domain over which α is calculated is shaded in green in figure 6. The model yields α = 2.69 on the beach, which is lower than in the coastal and  (2014) is used to determine the power law size range: for most scenarios, α is calculated over the entire size range except for scenario 3, which is shown using the green shading.
ocean cells (both α = 3.37). A scenario where the fragmentation rate is based on polypropylene instead of polyethylene is presented (red lines). Fragmentation breaks down the particles more quickly: a monotonic relation between particle size and abundance is observed, with α = 3.03. Finally, a scenario is presented (purple lines) where the input of plastic waste into the Mediterranean is 100,000 tonnes per year (Liubartseva et al., 2018;Jambeck et al., 2015), instead of the aforementioned 2,500 tonnes per year (Kaandorp et al. (2020)). The magnitude of the sinks needs to be much larger now to attain a mass balance based on 2,000 tonnes of floating plastics (Cózar et al., 2015). Fragmentation has little time to break down the particles, resulting in relatively few fragments per parent object. In figure 7, we compare PSDs resulting from the box model with observed ones. In the model results we include both size dependent ocean transport and resuspension. Fragmentation parameters are set to λ = 2 · 10 −4 f week −1 , and D N = 2.5, resulting in good agreement with the observed PSDs. The effect of vertical turbulent mixing of fragments using the model from Poulain et al. (2019) is shown as well (calm, U 10 ≈ 4 m/s, and above average, U 10 ≈ 7 m/s conditions based on the 30% and 70% quantile of Mediterranean sea weather conditions (Hersbach et al., 2020), see supplementary material section S3), assuming a submerged net depth of 25 centimeters. Figure 7a presents the resulting NSDs. Under calm wind and wave conditions there is good agreement between the modelled and observed NSDs. Vertical mixing causes the NSD to deviate from the power law around <3 milimeters, similar to the NSDs measured by Cózar et al. (2015) and Ruiz-Orejón et al. (2018). Many of the smaller fragments are expected to be mixed below the net depth, resulting in measuring only a fraction of small fragments. This, combined with a size detection limit effect where elongated particles escape from meshes smaller than their maximum length (Abeynayaka et al., 2020), could explain a part of the underabundance of sub-milimeter fragments in observations. Measurement campaigns with much smaller size-detection limits than the standard neuston nets (see e.g. Kooi and Koelmans (2019)) show increasing abundances for sub-milimeter fragments. It is therefore unlikely that the underabundance of sub-milimeter fragments is explained by an increased loss of these particles, suggested in some studies (Cózar et al., 2014;Pedrotti et al., 2016).
Resulting power law slopes α are presented in table 1, as well as the estimated power law size range. The model predicts a slightly larger α in the ocean compared to the coast. A similar difference is seen between the chosen sets of measurements near the coast (Ruiz-Orejón et al., 2018) and further away from the coast (Cózar et al., 2015), although this difference is not significant. Including vertical mixing has a strong effect on the estimated α.
Few PSD measurements are available for beaches. Two examples are shown in figure 7a: one from the Mediterranean (Constant et al., 2019), and one for which both the NSD and MSD were available (Fok et al., 2017). The measurements on beaches have much lower values for α compared to measurements in the water. This is also captured by the model, meaning size-selective resuspension indeed seems to play an important role. At the beach the modelled α values are higher than the measured ones, which might indicate that sizeselective beaching should be taken in account as well. Figure 7b presents the modelled MSDs. Vertical mixing has a large influence on the measured mass for small particle sizes: even under calm conditions, the measured mass for particles of 0.1 mm is almost three orders of magnitude lower than without mixing. Unfortunately, there is very limited observational data reporting MSDs, so the comparison to data is more limited than for the NSDs in figure 7a. On beaches, the model matches the set of measurement well, but more data are necessary to further verify this. Large fragments are expected to dominate in terms of mass on beaches. In the water, α seems to be approximately zero on average. This would mean that the mass contribution would scale roughly quadratically for an increasing size class k, i.e. large fragments also dominate in terms of mass here. The environmental box model used to model the PSDs is a useful tool for future mass balance studies. The steady-state with the model settings used for figure 7, gives that about 98% of the mass in the system is on the beach, about  This large fraction of plastics stranding is in good agreement with previous mass balance estimates (Lebreton et al., 2019). It should be noted that other environmental regions, like the ocean floor, are not included in these numbers as these are part of the sinks in the box model (P S ), which continuously take up more mass over time. Secondary microplastics generation can be estimated: for the same model settings, about 6.5 · 10 −5 % of the macroplastic (> 5 mm) mass breaks down into microplastics per week, about 2.0 · 10 −6 % of microplastics become smaller than 0.1 mm. This is orders of magnitude smaller than the estimated sinks, taking up about 5.0 · 10 −2 % of the plastic mass per week. Longevity of plastics can be estimated: taking a sudden stop of new plastics, it would take about 176 years for 99% of the plastic mass to disappear from the surface water and beaches. This is a much longer time scale than given by the conceptual model in , where for a similar stop of new plastics, almost all plastic mass was removed from the ocean surface layer within three years.

Conclusions and discussion
Modelling size distributions of marine plastics under the influence of continuous cascading fragmentationrtical mixing in the water column has a strong impact on measured PSDs. Understanding the nature of PSDs and how they differ in environmental regions can help close the marine plastic mass budget. It can be estimated which particle sizes contribute to most of the mass, as well as how persistent these are. The environmental box model is highly idealized using e.g. the assumption that the removal rate is the same in all compartments. Combined with more data and complex models (Kaandorp et al., 2020), this type of box model can be very valuable because different input and fragmentation scenarios can quickly be evaluated to assess whether these lead to realistic PSDs. The environmental box model was calibrated to Mediterranean sea data as much as possible. However, resuspension probabilities from Hinata et al. (2017) were estimated on beaches in Japan, and in Song et al. (2017) the experimental weathering was related to environmental time in South Korea. Further studies are necessary on how the parameters might vary globally.
The fragmentation model introduced in Cózar et al. (2014) focused mainly on spatial dimensionality: α = 3 in the NSD was related to three-dimensional fragmentation, i.e. a cube splitting into 8 smaller cubes. Care should be taken in future studies that when working with logarithmic binning, the normalized NSD (units: n mm −1 ) slope decreases by one compared to the discrete NSD (units: n), see figure 3. This was overlooked in Cózar et al. (2014): α = 3 would correspond to two-dimensional fragmentation with their model, see the supplementary material (section S5) for further explanation. Normalization is also important to take in account when describing plastic particle size in terms of a probability density function (Kooi and Koelmans, 2019), specifying the probability per unit length (units: mm −1 ). Finally, estimating α is not trivial: fitting straight lines on log-log transformed data induces large biases, maximum likelihood approaches are more suitable, see e.g. Newman (2005) and Virkar and Clauset (2014). The cascading fragmentation model by Charalambous (2015) used in this work, is able to model many small chips breaking away from a parent object, and showed good correspondence with experimental data from Song et al. (2017) (figure 5). In the cascading fragmentation model, the fraction of small particles keeps increasing over time, increasing α. We hypothesise that the value of α in the environment remains limited, due to the fact that smaller particles tend to be older, and hence are more likely to be removed by environmental sinks.
Future studies should further validate the cascading fragmentation model. We assumed that f and time are directly proportional here, which does not necessarily needs to be the case. In e.g. Charalambous (2015), it was shown that grinding can become less efficient as particles become smaller, which might lead to a logarithmic relation instead. Fragmentation of plastics is interesting in this sense: oxidation of plastics causes a brittle surface layer with microcracks (Andrady, 2011), perhaps increasing the fragmentation rate below a certain length scale, dependent on how far UV radiation penetrates the polymer.
MSD data is currently lacking: more data would help us to verify whether the larger size classes indeed make up most of the environmental plastic mass. More PSD data from beaches would allow for better constraining residence times of plastic particles on beaches and in coastal waters, and more data from marine sediment might give insight in the role of size-selective sinking, induced by e.g. biofouling .