A model for the size distribution of marine microplastics: A statistical mechanics approach

The size distribution of marine microplastics provides a fundamental data source for understanding the dispersal, break down, and biotic impacts of the microplastics in the ocean. The observed size distribution at the sea surface generally shows, from large to small sizes, a gradual increase followed by a rapid decrease. This decrease has led to the hypothesis that the smallest fragments are selectively removed by sinking or biological uptake. Here we propose a new model of size distribution, focusing on the fragmentation of marine plastics. The model is inspired by ideas from statistical mechanics. In this model, the original large plastic piece is broken into smaller pieces once by the application of “energy” or work by waves or other processes, under two assumptions, one that fragmentation into smaller pieces requires larger energy and the other that the occurrence probability of the “energy” exponentially decreases toward larger energy values. Our formula well reproduces observed size distributions over wide size ranges from micro- to mesoplastics. According to this model, the smallest fragments are fewer because large “energy” required to produce such small fragments occurs more rarely.


Introduction
A large fraction of the estimated billion tonnes of plastic waste that goes into the ocean [1] is found in a fragmented form, "microplastics", with sizes of less than 5 mm [2] through photodegradation and weathering processes [3][4][5].Those microplastics spread globally [6][7][8][9], potentially acting as a transport vector of chemical pollutants [10][11][12] and causing physical and chemical damages on marine organism [13][14][15].Recent drift simulations of microplastics calibrated against observed abundance of microplastics have produced global or semi-global maps of estimated microplastics abundance and concentration near the sea surface [6][7][8].These simulations will be further used to assess the biological impacts of microplastics.
Such simulations generally assume the size distribution of microplastics and their results would depend on the assumption because sedimentation and biological uptake can depend on size [14,[16][17][18].Interestingly, the number of small pieces rapidly decreases toward smaller sizes of O(100 μm) in most observations of plastic fragments at the sea surface [19][20][21][22].This feature is puzzling because the number of plastic pieces is expected to increase toward smaller sizes if the pieces keep broken down into smaller and smaller pieces (progressive fragmentation).For example, Co ´zar et al. [20] indicates that a type of progressive fragmentation leads to a cube law toward smaller sizes.The observed decrease at smaller sizes has generally been hypothesized to be due to selective sinking to depths, sampling error, or to selective ingestion by marine organisms [19][20][21][23][24][25][26][27].
The fracture mechanisms of marine plastics, however, are not well known.The power law, which results from scale-invariant fracture processes, is often invoked to explain observed size distributions of plastics.It tends to fit well observed size distributions at larger sizes [20].The power law can be derived, for example, from collision cascade among objects, as often applied to the fragmentation of asteroids [28][29][30], which does not include any decrease toward smaller sizes unless or until other processes than fragmentation start to dominate.Somewhat different fragmentation processes lead to a log-normal distribution [31][32][33], which is premised on the size reduction rate following the Normal distribution and has been applied to the fragmentation of sand grains, mastication, and the fracturing of thin glass rods [34][35][36], but not to oceanic microplastics to our knowledge.The log-normal distribution has a peak skewed toward smaller sizes and is similar to observed size distributions of fragments at smaller sizes but tends to deviate from observed distribution at larger sizes, where the power law tends to fit better [20].To the best of our knowledge, only the Weibull distribution [37,38] is similar to observed size distributions across the microplastics and mesoplastics ranges (fragments larger than 5 mm), at least in Co ´zar et al.'s study [20].The Weibull distribution as applied to the size distributions of various types of particles is largely empirical but with an interpretation as resulting from the branching tree of cracks [38].
The microplastics are thought to result from fragmentation by mechanical forces associated with waves or winds [3][4][5], but linkage between the microplastics and waves or winds has not been formulated in a physical manner.In this paper, assuming that "energy" needed to break down the plastic pieces acts as a constraint, we propose a new approach to modeling the fragmentation of marine plastics in relation to potential mechanical forces exerted on the plastic pieces borrowing ideas from statistical mechanics.As we shall see, under the assumption that the probability of the "energy" obeys an exponential distribution, the model reproduces a size distribution having a log-normal-like shape toward the small-size end and a power-law shape toward the large-size end.Both log-normal and power-law distributions have been commonly identified in observations from the surface to the mesopelagic levels in the ocean and on a beach [8, 19-22, 24, 39, 40].We will further discuss a potential explanation of the observed size distribution under the presence of deep intrusion or sampling error.

Fracture model
We assume that plastic waste first becomes fragile on beaches because of ultraviolet light and other weathering factors and then is broken up on beaches into microplastic pieces by the action of waves, winds, and other forces before being washed off into the ocean.Accordingly, the amount (mass) of microplastics which is produced on beaches will be determined by these processes.The fraction of microplastics that goes into the ocean may depend on waves and other weather factors.The microplastics will then be diluted in the ocean by mixing due to turbulence.These are the microplastics which observations sample.In the following we build a simple, idealized statistical model that determines the size distribution of the microplastics produced on the beaches.The model considers the fragmentation of given plastic mass and is naturally mass-conserving.
We offer the following simple physical model only as a representation of complicated fragmentation processes such as one-time crush and slow weathering.We do not claim that the following are exactly the underlying processes that lead to the observed size distribution.This model is originally intended to explain microplastics larger than O(100 μm).It will be later extended for finer microplastics (see the Discussion section and S4 Appendix in S1 File).In this section, we only show an outline of the derivation.Details are found in "Details of fracture model" of the Materials and Methods section.
The model assumes that the original plastic piece is a square plate with a size of L × L and a uniform thickness of Δh.This plate is broken into n × n equal-sized cells (Fig 1a).In this case, the size of each broken piece is λ = L/n.The number of the fragments of this particular size is given by n 2 multiplied by the number of the original plates (Fig 1a).
The fragmentation is assumed to be caused by "crush energy" that represents any mechanical forces such as winds and waves.In general, the energy required to fracture a brittle material can be evaluated by the surface energy (defined per unit area) multiplied by the area of the newly created surfaces (Griffith's theory [41,42]) (Fig 1b).Exploiting this fact, we define the crush energy here as the total surface energy of the broken pieces for each n.Highly plastic or nonlinear fracture processes involve other types of energy than surface energy [43][44][45], but we assume that breakage happens only after sufficient degradation and therefore that the plastic material is fully brittle [4,5], in which case, only surface energy needs to be considered.For a uniform surface energy density over the new surfaces, the total surface energy of broken pieces is proportional to the total length of the contact boundaries (Red lines in the right panel of  proportional to the total length of all contact boundaries.This means that to produce smaller pieces requires larger energy, which ultimately limits the number of small fragments as shall be seen below.This is the most significant new element of our model. Plastic pieces are fragmented by the action of waves, winds, or sand under various conditions, for example, on a hard reef or on soft sand, and therefore energy exerted on the plastic pieces is considered random.Unless each of these processes is modeled in detail, however, it is impossible to calculate the probability distribution of the energy from first principles.Nor is it possible to infer the probability distribution from field measurements at present.To make progress, we assume that an exponential distribution, p(ε) / e −ε/γ , governs the occurrence probability of crush energy ε (Fig 1c).This distribution is called "Boltzmann distribution" in the field of the statistical mechanics, which is often assumed by default when details of the underlying stochastic processes are not known and the distribution is known to be applicable to a wide variety of situations [46][47][48][49].This generality comes from the fact that it is the distribution that maximizes the entropy under fixed total energy [50]; that is, it is the most probable energy distribution under fixed total energy.We note that the long-term statistics of the coastal wave height generally shows that a larger wave height occurs with a smaller frequency [51][52][53], qualitatively consistent with the exponential distribution.
According to this probability distribution, a crush event with a large energy value is less frequent, consistent with common expectation.The factor γ may be regarded as a representative energy level of the environment, which represents the aggregate impacts of the combination of weather conditions (winds, waves, etc.) and the background conditions (hard or soft surfaces, etc.).
The number of the fragments of a particular size is n 2 j for a single crush event (Fig 1a).Since crush events occur randomly, the size spectrum is the expected value of n 2 j as a funcrtion of λ.Calculating the expected value with use of n = L/λ, we arrive at the size spectrum (see "Details of fracture model" of the Materials and methods section) where A is an arbitrary positive constant, which is adjusted by the amount of the total mass of plastics to be fractured.This size distribution has a shape skewed to smaller sizes similarly to the observed size distribution of the microplastics, as shown in Fig 2 .For a fixed A, the size distribution is controlled by γ and b.The factor b depends on the plastic material and the size of the original plate.The increase of γ and the decrease of b, however, have the same effect on controlling the size distribution and hence we introduce a new parameter γ � � γ/b.This parameter measures the strength of the mean environmental energy against the strength of the plastic material and its inverse provides a characteristic size such that the peak size, λ p , is given by λ p ' 0.255γ �−1 .As this parameter increases, thus, the size at which the maximum of the size distribution occurs decreases like γ �−1 and the corresponding maximum value increases like γ �4 (Fig 2a).Furthermore, in the large size limit, i.e., λγ � � 1, this size distribution asymptotes to S(λ)dλ * Aγ � /λ 3 dλ (Fig 2b), consistent with the cube power law observed in the mesoplastic range [20].
The total abundance and mass of the fragments for λ < Λ (�L) are obtained, respectively, as where σ � 2.404 (S1 Appendix in S1 File), ρ is the plastic density, and Λ is the maximum size of plastic fragments.Both of the above approximations are for γ � Λ � 1, that is, Λ being sufficiently larger than the characteristic size λ p .See S1 Appendix in S1 File for details.
We fit our model ( 1) to an observed size spectrum by adjusting A and γ � .Observed size distributions are presented in various ways: as the number per unit volume of sampled sea water, as the number per unit surface area of the ocean, as the raw number of collected fragments, etc.Some studies normalize the spectrum so that R S(λ)dλ = 1 [24,25,39].However, note that the size spectrum (1) can take any of these forms by adjusting A and different representations of the same spectra do not affect our analysis because the differences are absorbed into A.
We begin by comparing the present theory with the observed size distribution obtained by Isobe et al [21].This size distribution is based on the collection of plastic fragments sampled around Japan (Fig 3a ); this is the largest collection with the highest size resolution to date.For a precise comparison, we have converted the original size histogram into spectral density.Details are shown in "Observed data in Isobe et al" of the Materials and methods section.

Application to observed data
Theoretical curve is fitted to the observed size spectrum by adjusting the parameters A and γ � by a least-squares method over λ < 5 mm.The theoretical curve fits the observed size distribution well over the whole microplastic range (Fig 3b) with the relative error, Err = 7%.(The relative error is defined by the norm of difference between the theoretical and observed size distributions divided by the norm of the observed size distribution.)For comparison, we also plot a lognormal distribution, aLN(μ, σ 2 ), with an amplification or normalization factor a, where the parameters a, μ, and σ are determined by the same least-squares method.The lognormal distribution is not much different from our model in this size range.
Our model also agrees well with the observed size distributions in the Western-Pacific transoceanic section (Err = 22%) and the Seto Inland Sea (Err = 21%) (Fig 3c and 3d) with somewhat larger error than around Japan.The observed spectra are not as smooth as that around Japan, suggesting that the samples are not sufficient to give smooth distribution and this may be the reason for the larger error.Also, the peak is located at a smaller size in the Seto Inland Sea than in the other two regions; this shift is reflected in a larger value of γ � .The lognormal distribution is somewhat more different from our model in these two regions (Fig 3c and 3d) than in the area surrounding Japan (Fig 3b).The empirical fitting curve used in Isobe et al [8] has larger error, despite having the same number of parameters as ours (Fig 3c).
We next compare our model with another set of observations summarized by Co ´zar et al [20] (Fig 4a) based on samples collected in the accumulation zones ("garbage patches") over the globe, where microplastics originating from coastal areas tend to converge by ocean currents [55,56].Co ´zar et al's method of sampling plastic fragments is similar to Isobe et al's [8,19,21], but their size distributions include more size bins in the mesoplastic range than Isobe  and this discrepancy exists also in the case without the South Atlantic data (not shown).Despite having one more adjustable parameter, the log-normal distribution fit the data less well than our model.
The fitting of our model to observations gives a geographical distribution of γ � .The γ � value for Co ´zar et al's observation is largest in the Pacific Oceans, smallest in the North Atlantic, and inbetween elsewhere (Fig 4c).Given that the plastics are likely to be fragmented on beaches possibly by ocean waves [4], the value of γ � would represent wave energy on beaches where plastic waste is littered or comes ashore before fragmented there and washed away into the ocean as micro-or mesoplastics.A scenario-based numerical experiment (see "Wave energy in source region" of the Materials and methods section) suggests that the coastal area located along the western boundary of each basin is likely the major source region for the plastics in each accumulation zone, except that a large fraction of the plastics in the South Indian Ocean comes from South East Asia [55].Compared with a map of categorized wave energy levels along the world coastlines [57], the wave energy level in the potential source regions for a basin (Table 1) appears correlated with the γ � value for the accumulation zone of the basin in the emission scenario considering impervious area of land (Fig 4c; the results for the other scenario are found in Table 1 and S1 Fig in S1 File).The only exception is the North Pacific: the North Pacific accumulation zone has a large γ � value in spite of the relatively small wave energy level in the South China Sea, likely the major source region for the North Pacific.In fact, a large fraction of the plastic emission from the East China Sea is likely to come ashore on the southern coast of Japan [55], where the wave energy level is large (Table 1).If the majority of

Sources
Wave plastic waste is fragmented there, therefore, the large γ � value for the North Pacific may be explained (Fig 4c).The interpretation of γ � for Isobe et al's observations is less straightforward because the source region is less clear.The relatively low value of γ � for the area surrounding Japan (Fig 3b and S1 Table in S1 File) may be due to the low wave energy in the East China Sea, which may be the source region of the microplastics [55].In contrast, there are clearly multiple source regions contributing to the Western Pacific transoceanic section and it is possible that the distribution is a superposition of distributions with different γ � values (Fig B in S3 Appendix in S1 File).It is interesting that the Seto Inland Sea has a very large γ � (Fig 3d and S1 Table in S1 File), which may be due to some local conditions or to the conditions of some remote locations where the microplastics originate.It is indeed possible that some microplastics in this region originate from the Philippine Sea, where wave energy is large [57], as previous studies indicate that some waters from the Philippine Sea are transported into the Seto Inland Sea through the western boundary current [58][59][60].

Discussion
In summary, we derive a theoretical size distribution of micro-and meso-plastics using a statistical mechanical approach.It assumes that larger "energy" is required to break down the original plastic piece into smaller fragments and that this energy follows the Boltzmann distribution.This model well reproduces observed size distributions from the micro-to meso-plastics.In particular, it naturally explains the rapid decrease toward smaller sizes without invoking a removal of smaller fragments.This model is highly idealized, and extending this model for more realistic fragmentation processes that may be involved is a future study.
As pointed out above, Co ´zar et al's South Atlantic size distribution deviates from our theoretical curve more than the other size distributions in the same dataset.S3 Appendix in S1 File considers theoretical size distributions that would result if the sample is a mixture of plastic pieces originating from various source regions with different values of γ � , assuming that each source contributes the same number of plastic fragments for simplicity.The resultant superposition deviates from the single-source size distribution in a similar way that Co ´zar et al's South Atlantic size distribution deviates from the single-source theoretical distribution (Fig C in S3 Appendix in S1 File).Interestingly, the value of 1/γ � that best fits the mixture is close to the average of 1/γ � values of the origins in this simple example (See "Size distribution" in S3 Appendix in S1 File).
Similarly, a mixture of plastic fragments produced in various conditions may need to be considered for interpreting the observed total plastic mass.Isobe et al [8] show that the total mass of the plastics for sizes below 5 mm is about 0.3 mg m −3 for their samples in the Western-Pacific transoceanic surveys assuming that ρ = 1000 kg m −3 .To obtain this value from our M formula (2) with the optimal values of A and γ � for their dataset and with Λ = 5 mm, we arrive at a value Δh � 1 mm.Obviously this value of Δh is too large for the fragmentation model presented in this study because fragments smaller than Δh cannot be produced by the two-dimensional fragmentation of a plate with a thickness of Δh.This problem may be resolved if the value of Δh obtained like this is in fact an average of the various thicknesses of the original plastic plates (see "Total mass" in S3 Appendix in S1 File).
We have attempted to explain the observed size distributions at the surface by fragmentation alone, but in reality, a large fraction of the fragments are likely to be missing from observations.For example, deep intrusion is thought to be necessary to explain the amount of floating plastics [8,61], which is far less than estimates of terrestrial plastic emission [6,20].Biofouling, accumulation of organic matter on the surfaces of plastic fragments, is probably the primary mechanism for microplastics to lose buoyancy and sink to depths.This mechanism is more effective for smaller microplastics and therefore can affect the size distribution of microplastics at the sea surface [62][63][64].Sampling error can also affect observed size distributions.Some studies report that the rapid decrease toward the small size common to most size distributions at the surface may be an artifact attributable to the method of collection or size detection [24][25][26][27].If so many smaller microplastics are actually missing from observations, either by sinking or by escaping the nets, as to significantly alter the observed size spectra, the γ � values we have obtained by fitting are obviously underestimated, supposing, of course, that our theory is correct.
Further, secondary fragmentation may produce finer fragments.Since the smaller microplastics are closer to three-dimensional shapes than to a plate, their fracture would be a threedimensional process.The 3-dimensional version of our model (Fig. D in S4 Appendix in S1 File) can be similarly constructed to the 2-dimensional version (Fig 1a).The result is such that the model spectrum (1) gets one more factor of λ −1 and γ � = γ/(3L 03 ϕ), where L 0 is the size of the original microplastic cube.The peak size is at λ p = 0.201γ �−1 .See S4 Appendix in S1 File for a complete derivation.The theoretical curves of the 3-dimensional model are found to be in good agreement with the observed size spectra of the fine microplastics with scales of O (10 μm) collected near the sea surface [24,39], in the upper ocean from the surface to the mesopelagic layer [40], and on a beach [22] (Fig. E in S4 Appendix in S1 File).Those theoretical curves are characterized by γ � values of 1.68-14.01mm −1 , which are 10-100 times larger than the two-dimensional γ � values of 0.24-0.39mm −1 for Isobe et al.'s [8,19,21] and Co ´zar et al.'s [20] (Figs 3 and 4).Note that, however, the energy requirement for the fragmentation is given by γ.As an order-of-magnitude calculation, the value of γ = 3L 03 ϕγ � from the threedimensional γ � value is O( 103 ) smaller than the value of γ = 2L 2 Δhϕγ � from the two-dimensional γ � value under the reasonable assumptions that L = 100 mm and Δh = L 0 = 1 mm and that the surface energy density ϕ is the same.This small γ value comes from the smallness of the initial size, L 0 , which we assume that the three-dimensional fragmentation starts from.Plastics with smaller initial sizes can be fragmented by smaller energy since the surface energy decreases with the initial size (see Fig 1b  Accordingly, the fate of the marine plastics may be as follows.Original plastic pieces deposited at a beach are fragmented down to millimeter size by weather phenomena such as waves.The two-dimensional version of our model intends to predict the spectrum of these fragments.Those microplastics are further fragmented into finer pieces by slow weathering, grinding in sand, or other processes.The three-dimensional version of our model may also explain the resultant size spectrum of the finer fragments.Both types of microplastics are washed off into the ocean.Given that a smaller particles are susceptible to the vertical mixing [19,23,24,65], those finer microplastics may be spread over the mixed layer or deeper, while the larger ones may tend to stay at the sea surface.This scenario based on our fracture models has no contradiction to the observed evidence on the spatially varying size distribution of the microplastics so far. The reduction in the particle count at smallest sizes would be a physically reasonable feature when available fracture energy is limited.Indeed, the general grinding needs larger energy to create smaller particles [66,67], and the size distribution resulting from a progressive grinding process has this feature [68].The power law, which the literature often invokes to explain size distributions, is usually ascribed to some progressive fragmentation processes.Such processes may also result in a decrease at small sizes if a limitation of energy is introduced as it is to our model.When the observed spectrum drops at smaller sizes, the lognormal distribution is also often invoked and contrasted to the power law in the literature on experiments for plastic fragmentation [69,70].Quantitatively, the lognormal distribution does not much differ from our model spectrum (Figs 3 and 4) although it does not agree with the power law in the larger-size range (Fig 4d).Distributions which were compared to the lognormal distribution in the literature may be explained by similar mechanisms to our model.The present model is premised on the fragmentation of a plate into equal-sized cells with a regular square shape, which provides the inverse proportionality in the crush energy to the size (See section).However, this square lattice is merely a simple device to relate the total surface energy (the work required to create new surfaces) with the size and the number of fragments, but this regularity is not necessary.Indeed, the same relationship between the crush energy and the size holds for irregular cells given by a Voronoi tessellation [71] (S5 Appendix in S1 File), which is often used for modeling pervasive fracture such that an impact causes propagation of multiple cracks and collapse [72][73][74][75][76].The pervasive fracture caused by an instantaneous impact generally can yield the fragments with various sizes.Our modeled size spectrum may be regarded as an ensemble of the fragment counts for each size for multiple pervasive fracture events with a nearly uniform impact.
The scale invariant approach treating fractality in fracture is often used to describe the fragmentation process of materials into isolated pieces and can explain power law size distributions in fragment counts observed in the breakage of the materials such as concrete, glass, and ice in nature [77][78][79].In these analyses, consistent with Griffith's theory of crack extension, the energy consumption to create fragments are represented by surface energy, which is proportional to the surface area of the fractures [80][81][82][83][84].Our model shares this property.An essential difference between the two approaches is that our model limits the available energy by introducing occurrence probability of crush energy.As a result, our distribution takes a qualitatively similar form to the log-normal distribution in the smallest-size region and reduces to a power law in the large-size limit (See Fracture model section).The latter limit is given by λγ � � 1, which is equivalent to the large γ � limit; that is, when γ � is large, the power law region extends leftward (toward smaller sizes in the size spectrum).Unlike in fractal processes, the exponent in the power law size distribution is an integer number in our model, which is a consequence of crush energy (total surface energy from a single plastic plate) being inversely proportional to the size of the fragments.In the above argument, we have shown that this proportionality holds even for the irregular fragment cells following a Voronoi tessellation but we have not permitted surface fractality.Further generalization of our model is needed to cover a wide range of fractal processes.
The present fracture model may open paths toward developing sophisticated numerical simulations for predicting the production and dispersion of the microplastics.There have been numerical simulations to evaluate the spreading microplastics in the ocean [6][7][8][9]85].In such simulations, virtual parcels representing a mass of plastic pieces are released at source regions and advected by ocean currents.The amount of plastics released has to be either assumed or calibrated so that the resultant mass distribution matches observations.The size distribution, if that information is necessary, is usually assumed in an ad-hoc manner or calibrated against observations.Given that the microplastics are likely to originate from the weathering of plastic litter on beaches [3][4][5], their size distribution would depend on weather and wave conditions that would be different for each beach.Our size distribution model may be used to estimate the initial size distribution of plastic fragments by parameterizing γ � at the beaches as a function of weather conditions such as winds and waves.

Details of fracture model
In the Fracture model section, we only outlined the derivation of our model spectrum and focused on its interpretation; here we show full details of the derivation.Our model is derived in analogy with black body radiation; see S2 Appendix in S1 File for the analogy.
In this model, we consider that a square plastic plate with a size of L × L and a uniform thickness of Δh is broken into n × n equal-sized cells (Fig 1a).The size of each broken piece is given by λ = L/n; we also define an associated "wavenumber" ν � λ −1 = n/L for convenience.
In general, the minimal energy required to fracture a lump of solid is given by the surface energy [86], which is proportional to the area of the newly created surface.This area in the present model is proportional to the total length of the contact boundaries (Fig 1a), which is related to n, like 2L(n − 1) � 2Ln.With the aid of ν = n/L, the crush energy for a plate can therefore be written as bν with a constant b � 2L 2 Δhϕ, where ϕ is a uniform surface energy density, and hence the crush energy for j plates can be written as ε � jbn: ð3Þ We assume that the occurrence probability of crush energy is governed by the Boltzmann distribution characterized by γ: which states that a crush event with a large energy value is less frequent, consistent with common expectation.In statistical mechanics, γ is temperature in heat bath times the Boltzmann constant (Fig A in S2 Appendix in S1 File).This is the analogy to the "environment" for microplasitics.From (3), the expected value of the crush energy as a function of ν, denoted by hεi ν , is and from (4), hji n ¼ which is the expected value of the number of fractured plates for each ν.In deriving (5), we used the geometric series such that S 1 j¼0 ar j ¼ a=ð1 À rÞ, which holds for any a and any r 6 ¼ 1.This formula is called the Bose distribution [87].Since the number of the fragments of a particular size is n 2 j (Fig 1a ) and n = Lν by definition, the expected number of fragments is P(ν) / n 2 hji ν / ν 2 hji ν and therefore, or converting from ν to λ, we arrive at the size spectrum we already presented in the Fracture model section where A is an arbitrary positive constant.Eqs 6 or (1) is formally analogous to the Planck's formula of blackbody radiation (See S2 Appendix in S1 File).

Observed data in Isobe et al
The observed size distributions in Isobe et al are from surveys around Japan [21], in a western Pacific transoceanic section [8], and in Seto Inland Sea [19] (Fig 3a).The first set of samples is from 56 stations around Japan during the period of July 17 through September 2, 2014 [21].
The second is a transoceanic survey at 38 stations across a meridional transect from the Southern Ocean to Japan during 2016 [8]; the stations in Southern Ocean and the other ones were occupied from January 30 to February 4 and from February 12 to March 2, respectively.The third dataset is collected at 15 stations in the Seto Inland Sea of Japan in May-September from 2010 to 2012 [19].For all surveys, neuston nets with mouth, length, and mesh sizes of 0.75 × 0.75 m 2 , 3 m, 0.35 mm, respectively, were used to sample small plastic fragments.The nets were towed near the sea surface around each station for 20 min at a constant speed of 2-3 knots (1-1.5 m/s).The numbers of fragments are counted for each size bin with a bin width of 0.1 mm for microplastics and 1 mm for mesoplastics (defined to be >5 mm) between 5 and 10 mm, except for the surveys in the Seto Inland Sea, in which the bin width becomes wider beyond the size of 4 mm.The size of a fragment is defined by the longest dimension of its irregular shape.The concentration of the fragments (pieces per unit volume of sea water) within each size bin were calculated by dividing the number of fragments by the water volume measured by the flow meter at each sampling station.This concentration binned according to size is the observational data used in the present study.To obtain the numbers, we have digitized the plots of Isobe et al's using WebPlotDigitizer version 4.3 (https://automeris.io/WebPlotDigitizer/).The original size distributions in the Seto Inland Sea are presented separately for four different areas [19], but they are averaged in the present study.
Fig 5a replots one of these size distributions as an example.As stated above, the bin width is not uniform: for this particualr data, Δλ = 0.1 mm for λ < 5 mm, Δλ = 1 mm for 5 mm < λ < 10 mm, and Δλ = 10 mm for λ > 10 mm.This is the reason that the concentration jumps up beyond λ > 5 mm.For comparison with theories, we introduce a "size spectral density" n i (Fig 5b) such that where Δλ i is the width of the i-th bin and N i is Isobe et al.'s value for the bin.By this definition, each rectangular area of the spectral plot is proportional to the number of plastic fragments within the bin.The spectral density is less sensitive to the bin widths, and in the limit of Δλ !0, it converges to a continuous size spectrum S(λ) such that R l b l a dl SðlÞ is the number of fragments between λ a and λ b .All of Isobe et al's size distributions are converted to size spectral densities in the present study.

Observed data in Co ´zar et al
The plastic samples summarized in Co ´zar et al [20] were collected in the accumulation zones [56]

Wave energy in source region
We use the Lebreton et al's plastic dispersal simulation result [55] to identify the source regions of the plastic fragments sampled in the Co ´zar et al's accumulation zones.Lebrenton et al's numerical simulation provides the dispersion of virtual particles representing a set of plastic fragments originating from land following the sea-surface currents reproduced by an oceanic circulation modelling system HYCOM/NCODA [88].The particles are released at coastal locations on the basis of two scenarios, the amount of released particles determined as a function of impervious surface area (ISA-based scenario) and of coastal population density (PD-based scenario), respectively.The former is intended to reflect contributions from major rivers and the latter from large cities.Contribution in percentage of each emission region to the amount of the plastics in the accumulation zones are summarized in Table 1 from Lebreton et al's supplementary data.For example, the plastic emission from Australia accounts for about 9.5% of the amount of plastics found in the South Pacific accumulation zone for the ISA scenario.
Next, we estimate wave energy level (no units) at each emission region from the 6-grade wave energy levels along the global coastlines provided by Fairley et al [57].Each of Lebreton et al's emission regions for each emission scenario consists of multiple emission points (their S2-S7 Figs).We select the point with the largest emission within the region, look up Fairley et al's figure to determine the wave energy level of the point, and assign the energy level to the region.If there are multiple points with similarly large energy levels, we assign the average energy level to the region.These energy levels are summarized in the "Energy Level" column of Table 1; the left and right subcolumns are for the ISA and PD senarios, respectively.Some of the energy level values are missing in the table because the selected emission points are not covered in Fairley et al's map.Finally, the expected wave energy level that the plastics found in the accumulation zone experienced for each scenario is calculated by the contributionweighted average ∑ i E i f i /∑ i f i , where E i and f i denote the wave energy level and the contribution of the i-th source region, respectively.For the North Pacific, we calculate the expected wave energy level without contribution from China, assuming that plastic wastes from China go to the North Pacific accumulation zone via Japan.The resultant expected wave energy level for each accumulation zone is shown in S1 Fig in S1 File.
Fig 1a) since the plate is assumed to have a uniform thickness.Because the total length of the contact boundaries increases linearly with n (Fig 1a, right), the required crush energy is

Fig 1 .
Fig 1.Schematic illustration of model configuration.a. Schematic representation of our fracture model.An idealized plastic plate with a size of L × L and a thickness of Δh is broken into n × n equal-sized pieces.The size of each broken piece is λ = L/n (left).The total length of the contact boundaries (red lines) is obviously l � 2(n − 1)L (right).We sum up energy needed to break j plates into λsized pieces and denote it as ε (middle).b.Schematic representation of breakage.Energy needed to break up the plate is proportional to the cross-sectional area of the contact boundary and hence to the length of the fracture.c.Occurrence probability of crush energy governed by the Boltzmann distribution, p(ε) = e −ε/γ /γ, which exponentially decreases for larger crush energy ε.The e-folding scale is γ and hence the probability of larger crush energy becomes larger as γ increases.https://doi.org/10.1371/journal.pone.0259781.g001

Fig 2 .
Fig 2. Theoretical size spectrum.a. size distributions expected from (1) for different values of γ � under a fixed A (= 1.0).The dashed curve connects the peaks of the size distributions.b. a log-log plot of the same distribution for γ � = 0.4 and A = 1.0 (solid curve).The dashed line denotes the power law Aγ � /λ 3 , which the size-distribution curve asymptotes to for large λ.https://doi.org/10.1371/journal.pone.0259781.g002 et al's (See "Observed data in Co ´" of the Materials and methods section).Our model generally reproduces well the observed size distributions although it has large relative error in the South Pacific Ocean (Fig 4b), where the least number of samples were collected in Co ´zar et al.'s study and the spectrum is the least smooth (barely visible in the plot), suggesting that the number of samples is not sufficient.The unweighted sum of the five size distributions is shown in Fig 4d in a log-log form.Importantly, the model reproduces the cube power law toward the mesoplastic range (Blue curve in Fig 4d) as found by Co ´zar et al.The theoretical curve is not particularly

Fig 3 .
Fig 3. Comparison of theoretical size spectrum with the observed size distributions of Isobe et al. [8, 19, 21].a. Schematic map of observation stations in Isobe et al. [8, 19, 21].See their original papers for the detailed locations.b-d.Theoretical and observed size spectral densities of microplasitcs in the area surrounding Japan, along a Western-Pacific transoceanic section, and in the Seto Inland Sea.The black bars display the observed size spectral density (left axis), and gray bars are the original histogram (right axis) taken from Isobe at al.In the range λ < 5 mm, the black and gray rectangles perfectly coincide.See "Observed data in Isobe et al" in the Method section for the detailed method of the conversion from the histogram to the spectrum.The blue and orange curves denote our model and the lognormal distribution, respectively.The parameter A is dimensionless (×10 −9 ) in this case.The green dashed curve in c denotes the empirical curve βλe −αλ of Isobe et al [8].The map in a is created using a Python package Cartopy [54] (version 0.18) and Adobe Illustrator CS6 (version 691; https://www.adobe.com).https://doi.org/10.1371/journal.pone.0259781.g003

Fig 4 .
Fig 4. Application of theoretical model to Co ´zar et al. [20]'s observational data and comparison between γ � and wave energy.a. Schematic map of observation stations in Co ´zar et al. [20].See the original paper for the detailed locations.b.Theoretical (curves) and observed (symbols) size spectral density of all samples for each accumulation zone.The observed data are same as those in their S6 Fig, except that the values are expressed as a spectral density (see "Observed data in Co ´zar et al" in the Materials and methods section).c.Expected wave energy levels (no units) at source regions (rectangles) and γ � (dots) for accumulation zones.The accumulation zones are denoted by abbreviations such as SI for South Indian Ocean.The red dashed rectangle shows the wave energy level for the North Pacific Ocean in the case where the contribution from China is removed.See "Wave energy in source region" in the Materials and methods section for details.d.Sum of the size spectral densities plotted in b over the basins for observation (black dots) and our model (blue curve).The orange curve and black dashed line denote lognormal distribution and a cube power law, respectively.The observed data and cube power law line are the same as those in Co ´zar et al's S10 Fig.The gray dots and curve are the same as the black dots and blue curve, respectively, except that the South Atlantic data is excluded for the sizes less than 30 mm (the digitizer does not recognize the small spectral values above 30 mm for the South Atlantic data in Co ´zar et al's figure).The scale of each axis on panels b and d follows that of Co ´zar et al's S6 and S10 Figs, respectively, except that the vertical scale on panel b is shown as spectral density.The map in a is created using the same softwares as in Fig 3a.https://doi.org/10.1371/journal.pone.0259781.g004

Table 1 .
Expected wave energy levels (based on Fairley et al's [57] map) at source regions and plastic emissions (percentage) that contribute to the five major oceanic accumulation zones based on the two terrestrial release scenarios (ISA scenario and PD scenario) from Lebreton et al. [55].The categorization of the source regions is from Lebreton et al's supplementary material.The values for the ISA and PD scenarios are indicated in the left and right sub-columns of each column, respectively.See the Methods section for the details.
in this text and Fig D in S4 Appendix in S1 File).

Fig 5 .
Fig 5. Size distribution expressed as a histogram (a) and as a spectral density (b) from Isobe et al's [21] observation around Japan.The histogram is a replot of Isobe et al's Fig 2. We have obtained the data by digitizing the original figure using WebPlotDigitizer version 4.3 (https://automeris.io/WebPlotDigitizer/).The size spectral density is indicated by black bars with its scale shown on the left axis of panel (b), and the gray bars on panel (b) are the same histogram with its scale on the right axis.The spectral density is plotted in such a way that the black bars exactly coincide with the gray ones for λ < 5 mm.In panel (b), sizes larger than 10 mm are omitted because the spectral values are almost zero there.https://doi.org/10.1371/journal.pone.0259781.g005 around the world (Fig 4a) from December 2010 to July 2011.Their method of sampling was the same as that of Isobe et al's studies, except for the different mesh size (0.2 mm at the minimum) and mouth area (1 × 0.5 m 2 ) of the neuston net and the different towing durations (10-15 min).The collected plastic fragments are classified into bins whose widths increase exponentially (like Δλ k = c10 0.1k mm, where k is the bin number) for 0.2 mm < λ < 100 mm.Co ´zar et al summed the number of plastic pieces over each accumulation zone without dividing the number by the volume of the sampled sea water (S6 Fig in their paper).Their data are digitized and converted to spectral densities in the same way as described above for Isobe et al's data.Sizes larger than 40 mm in the tail of the distribution have had to be omitted because the numbers are so small there that the data points are too close to the horizontal axis of the plots for the digitizer to resolve.This limitation does not significantly influence the curve fitting because the fitting result is most highly sensitive to the values of the size spectral density around the peak size.The present study has also digitized Co ´zar et al's S10b Fig to create a logarithmic plot of the sum of the size distributions for all the observed accumulation zones (Fig 4d).