A coupled Lagrangian-Eulerian model for microplastics as vectors of contaminants applied to the Mediterranean Sea

The pervasiveness of microplastics in global oceans is raising concern about their impacts on organisms. While quantifying their toxicity is still an open issue, sampling evidence has shown that rarely are marine microplastics found clean; rather, they are often contaminated by other types of chemical pollutants, some known to be harmful to biota and humans. To provide a first tool for assessing the role of microplastics as vectors of plastic-related organic pollutants (PROPs), we developed a data-informed simulation model that accounts for the intertwined dynamics of Lagrangian microplastic particles transported by surface currents and the Eulerian advection-diffusion of pollutants that partition on them through seawater-particle interaction. Focusing on the Mediterranean Sea and using simple, yet realistic forcings for the input of PROPs, our simulations highlight that microplastics can mediate the export of PROPs across different marine regions. Particle origin, in terms of both source type (either coastal, riverine, or fishing-derived) and geographical location, plays a major role in determining the amount of PROPs conveyed by microplastics during their journey at sea. We argue that quantitative numerical modelling approaches can be focal to shed some light on the vast spatiotemporal scales of microplastics-PROPs interaction, complementary to much-needed field investigation.


Introduction
Microplastics (plastic fragments smaller than 5 mm, Arthur et al 2008) are emerging as a new class of marine contaminants . Not only are microplastics composed of various types of polymers and chemical additives, as diversified as the products they originated from, but their chemical profiles may also be modified during their journey at sea. For instance, microplastics can leach the additives used in the production of their parent material (Koelmans et al 2014); they also have high sorptive capacities for hydrophobic environmental contaminants (Hartmann et al 2017) due to chemical affinity (O'Connor et al 2016). The mounting evidence about microplastics' relevance as vectors of contaminants is raising concern (Jiménez-Skrzypek et al 2021), because such microscale interactions may determine toxicity to organisms, up to impacting processes throughout the marine ecosystem (Galloway et al 2017). The negative impacts of microplastics on marine biota, such as blockage and/or reduced nutrition (Wright et al 2013), or oxidative damage (Galloway and Lewis 2016), up to translocation in tissues (see e.g. Mattsson et al (2017) may in fact be exacerbated by chemical risk. Microplastics could function as an additional exposure medium to environmental pollutants on top of other pathways, with potential synergistic toxic effects on organisms when combined with harmful chemicals-a point of concern still at the frontier of relevant research (Amelia et al 2021).
An important effort to provide quantitative insights about the pollutants carried by marine plastic litter is that of De Frond et al (2019), who combined field data on marine plastic litter with information about plastic manufacture to quantify the plasticmediated influx of additives to the global oceans and the amount of polychlorinated biphenyls (a typical pollutant with affinity to plastics) adsorbed on microplastics at four selected locations. Nonetheless, they did not track the consequent fate at sea of plastics, nor how and where pollutants were adsorbed on microplastic particles. Lagrangian particle-tracking models at oceanographic scales are widely used to describe microplastic transport and distribution on the sea surface of the global ocean (e.g. Lebreton et al 2012, among the pioneering ones) and of the Mediterranean Sea (where Mansui et al 2015 paved the way), complementing the insights obtained from field research (Hardesty et al 2017). Nonetheless, never have particle-seawater chemical interactions occurring during the journey at sea of plastic fragments been studied through a modelling framework. These chemical exchanges add a new dimension of complexity to the modelling of microplastic pollution, i.e. accounting for environmentally relevant changes caused by chemical-physical interactions to both microplastics and seawater superimposed to (a) the journey of each fragment and (b) the advectiondiffusion dynamics of plastic-related organic pollutants (PROPs, Guerrini et al 2021).
Here we propose a novel modelling framework for microplastic-related pollution combining Lagrangian tracking of microplastics and Eulerian advectiondiffusion of PROPs; these two sides of the model are coupled through microplastic-seawater chemical exchanges. Aim of our work is to describe the interlinked spatiotemporal dynamics of microplastics and PROPs within ocean the surface layer and apply here our methodology to the Mediterranean Sea. This semi-enclosed, highly urbanized basin is considered as one of the global ecoregions most affected by anthropogenic pressure (Halpern et al 2008); it has also been regarded as an accumulation area for plastic debris (Cózar et al 2015) and chemical pollutants (Berrojalbiz et al 2011). Within the Mediterranean, we identify areas where the exchange of PROPs operated by microplastics modulates seawater pollution. Microplastic-mediated pollution is also apportioned in relation to particle source-type (coastal, riverine, or fishing-related) and geographic area of origin. Additionally, we analyze particle trajectories to detect recurrent spatial patterns possibly driving microplastics-PROPs interactions at a basin scale.

Methods
We identified three separate Sub-Problems (SPs) to be solved when considering the coupled dynamics of microplastic particles and PROPs in marine environments (figure 1(a)): • SP1-particle transport by surface currents, • SP2-advection-diffusion of PROPs on the sea surface, and • SP3-mass transfer of PROPs between particles and seawater.
Numerical simulations of our coupled algorithm were designed using an operator-splitting approach (sensu Strang 1968), whereby each SP was solved over staggered time steps. More precisely, the temporal sequence of subproblems execution in our algorithm was defined as follows: (a) Eulerian advection-diffusion of PROPs over the oceanographic domain (SP2) is solved over a half time-step [0; ∆t 2 ]; (b) Lagrangian particle transport by surface currents (SP1) is simulated at a Mediterranean-wide scale over an entire time step [0; ∆t]; (c) Coherently with the chemical gradient between the concentration of PROPs on each particle and in the surrounding seawater at time t = ∆t, water-particle pollutant exchanges are solved at particle scale and concentration of PROPs in seawater is updated over [0; ∆t] (SP3); (d) SP2 is solved again for another half time step [ ∆t 2 ; ∆t], so that the numerical scheme is completed over a full time step.
This algorithm was then executed throughout the duration of the simulation.

SP1: Lagrangian particle tracking
We exploited the data-informed Lagrangian model presented in Guerrini et al (2021), briefly recalled in the following (further details in the supplementary material, SM available online at stacks.iop.org/ERL/ 17/024038/mmedia). As customary, particles were treated in SP1 as point-like entities moving coherently with surface marine currents. Drivers of particle movement were the velocity vectors obtained at particle coordinates by interpolating with a tri-linear scheme the zonal and meridional components of the daily currents from the surface layer of the Mediterranean Sea Physical Reanalysis (1/16 • × 1/16 • horizontal resolution and 72 layers of depth; (Simoncelli et al 2014)) available through Copernicus Marine Environment Monitoring Service. We did not overimpose to reanalyses any additional component of particle motion, as wind drag or wave entrainment, also because of the small size of our modelled particles. It is worth remarking that particles' motion do not depend on chemical exchanges; instead, the chemical state of seawater is influenced by that of particles. Thus, SP1 could be run first and with no reference to the other SPs.
The initial vertical position of each released particle was assigned randomly upon release within the uppermost layer of the reanalyses (1.5-4.6 meter deep), and was then kept constant throughout the simulation, without accounting for vertical velocities. The removal of particles from the sea surface is due to natural processes, here referred to as 'sinking' , caused by e.g. decrease in buoyancy due to the growth of biofilm, see e.g. Lobelle et al (2021) and/or aggregation in marine snow (Kvale et al 2020). It was parameterized by attributing to each particle a transport duration randomly extracted from an exponential distribution with an average of 50 days (Liubartseva et al 2018) and not exceeding 250 days. This resulted in a median particle sinking time of 34 days, a value coherent with the results by Fazey and Ryan (2016) and Lobelle et al (2021). Another mechanism acting on the balance of microplastics is beaching (Eunomia 2016, Kaandorp et al 2020). In our model, a particle was considered beached (and its simulation stopped) if it reached a position where both the zonal and meridional components of the velocity field were null. Resuspension of beached particles (e.g. because of rough sea or tidal motion) was disregarded.
New particles were released daily from significant sources of marine plastic pollution, of either landbased origin, as coastlines (Jambeck et al 2015) and rivers (Lebreton et al 2017, Boucher andBillard 2020), or offshore origin, namely fishing grounds (Ocean Conservancy 2011). We tracked around 100 million particles per year, assigned to each specific source by adopting the 50:30:20 coasts-to-rivers-to-ships ratio proposed by Liubartseva et al (2018). We further modulated in space and time the release of particles from each source point by adopting data-informed, source-specific proxies, as briefly detailed below (see also SM, section S1).
The mismanaged plastic waste (MPW) produced along the Mediterranean coastlines was selected as a proxy indicator of coastal microplastic input. Coastal sources thus emitted a number of particles proportional to the product between two factors: the country-specific estimate of per capita MPW generation rate (from Jambeck et al 2015) and the population inhabiting in a 5 km-wide strip along the coast (from CIESIN 2018). The proxy indicator of particle release from rivers accounted instead for MPW generated within drainage basins and for seasonal variations in surface water runoff, as the hydrological regime strongly affects the presence of debris in river waters , van Emmerik et al 2018. It was obtained by summing over each basin the gridded product of (a) monthly surface runoff (from the gridded 2015-2016 monthly FLDAS dataset, (McNally and NASA/GSFC/HSL 2018)), (b) inhabiting population (from CIESIN 2018), and (c) country-specific estimated per-capita MPW generation rates, using the same gridding as the CIESIN population data (the finest among the three datasets involved in the operation). We then averaged the indicator over 2015-2016 and selected as particle sources for our model the rivers ranking in the top-100, and modulated their particle releases according to the monthly values of their proxy indicator.
Microplastic input from fishing activities was modulated proportionally to the daily fishing effort at offshore locations from the Global Fishing Watch dataset (https://globalfishingwatch.org/; Kroodsma et al 2018) for the years 2015-2016. We removed data regarding fixed fishing gear (like longlines, pots, and traps), as we assume that the release of microplastics happens mostly during handling.

SP2: Eulerian advection-diffusion
To model the dispersion of PROPs within the surface layer of the Mediterranean Sea, we used a simple Eulerian approach based on the general advectiondiffusion equation (Fischer et al 1979): where C W represents concentration of PROPs in seawater, D is the diffusivity coefficient (table 1), v is the velocity vector, and S is a location-dependent input term (SM, section S2). To numerically solve  (1), we used the cell-centered finite volume method implemented in the PDE solver FiPy (www.ctcms.nist.gov/fipy/, (Guyer et al 2009)). For spatial coherence with SP1, we defined a control volume with the same grid size (1/16 • × 1/16 • ) and depth (3.1 m) used there. At present, data scarcity prevented us from informing the inputs of PROPs with field observations of hydrophobic pollutants at a whole-Mediterranean scale (Guerrini et al 2021). To parameterize the source term of equation (1), in SP2 we released a target PROP (phenanthrene) daily from the same sources as in SP1, following the proxies that guided the release of microplastics. Despite its simplicity, our Eulerian module proved to yield realistic concentrations of our target PROP on the sea surface (Guerrini et al 2021, and SM, section S2). Contrary to what remarked above for SP1, the chemical exchanges occurring at the particle scale (see SP3 below) do affect the concentration of PROPs in seawater, thus they interact with the advection-diffusion processes occurring at basin scale described in SP2. Therefore, SP2 and SP3, both informed by SP1, had to be executed at right timing (see the coupling algorithm above).

SP3: particle-seawater mass transfer
To quantify the mass transfer of pollutants between microplastics and seawater, particles transported in SP1 were qualified in SP3 with specific properties (like polymer type, prototypical size, and shape of microplastics; Rochman 2015) to (a) allow realistic microplastic-seawater partitioning of PROPs, and (b) more accurately describe transfer kinetics. As for (a), we assumed our particles to be 1 mm polyethylene spheres, and adopted a super-individual approach (sensu Scheffer et al 1995) to simulate the PROP-exchange capability of realistic concentrations of microplastics on the surface of the Mediterranean Sea (see SM, section S3), i.e. each particle represents a certain abundance of microplastics. Under our assumptions, the exchange of PROPs could be modelled by applying Fick's law to each particle, i.e.
where C i,j represents the concentration of PROPs on particle i within cell j, whose concentration of PROPs in seawater computed in SP2 is C W,j . Both the plasticwater partition coefficient K P−W and the desorption rate k des (table 1) depend on PROP-polymer pairing (SM, section S3.2). The concentration of PROPs on all particles i's are updated synchronously, i.e. at constant C i,j , an assumption made plausible by the widely different spatial scales at which SP2 and SP3 occur. Equation (2) was solved numerically for each particle using the Python package Scipy (https://docs.scipy.org/doc) with an explicit Runge-Kutta method of order 5(4). The concentration of PROPs in seawater was updated by summing up all the mass of PROPs exchanged (i.e. adsorbed or desorbed) between particles and seawater in cell j (SM, section S3).

Numerical experiments
By properly intertwining SP1, SP2, and SP3 at daily temporal resolution, we simulated two years (2015-2016) of coupled Lagrangian-Eulerian interactions between particles transported by Mediterranean surface currents and a target PROP undergoing advection-diffusion on the sea surface layer. To set credible, basin-wide initial conditions for the coupled simulation, SP1 and SP2 required a warm-up. As for SP1, we began the Lagrangian simulation starting from a sea surface clean of microplastic particles 250 days before 1st January 2015. This duration corresponded to the longest residence time of particles floating within the sea surface layer, after which the Lagrangian part of the model is effectively memoryless. All particles were assumed to be uncontaminated at release. Since the advection-diffusion of PROPs is not memoryless, we used a longer warm-up period (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) for SP2, starting with a clean Mediterranean Sea and performing daily releases of the target PROP from each source type and location (SM, section S2). During the whole coupled Lagrangian-Eulerian modelling simulation, we tracked trajectories, particle-seawater PROP fluxes, and daily PROP concentrations of more than 250 million microplastic particles. We then mapped simulation outputs, i.e. the time-averaged PROP concentrations in seawater and on particles, and the average direction (adsorption or desorption) and magnitude of particle-seawater exchanges with respect to particles (additional details can be found in SM, section S4). To further analyze the spatiotemporal patterns of microplastic-mediated PROP dispersal by apportioning particles (and the PROP conveyed by them) to their geographic origin, we grouped the maritime boundaries of each of the Mediterranean-facing countries (as in Marine Boundaries v11, Flanders Marine Institute 2019) into seven regions, as shown in figure 1(b). We then assessed whether highly PROP-polluted particles tend to accumulate and where. Highly polluted particles were identified as those exceeding a suitable ecotoxicological threshold of PROP concentration. To the authors' best knowledge, no environmental risk limits have been standardized yet for PROPs sorbed on microplastic particles; we thus relied on the maximum permissible concentration (MPC) for phenanthrene (our target PROP) in saltwater sediment proposed by Verbruggen and van Herwijnen (2011), corresponding to 780 ng PROP /g microplastic or higher. We then computed the temporally averaged concentration (in terms of no. particles km −2 ) of all particles whose pollution exceeded the MPC to map the hazard caused by particle-bound PROP in the Mediterranean Sea.

Results
Maps of the temporally averaged concentration (2015-2016 average) of the target PROP, as resulting from by our coupled Lagrangian-Eulerian model, are shown in figure 2. PROP concentration in surface seawater (figure 2(a)) follows a rather smooth, west-to-east increasing gradient. The most notable local exceptions to this trend (see figure 1(b) for reference) are the Gulf of Lyon (Northwestern region), the Albanian and Ionian Greek coasts (Adriatic region), and the southern Sicilian shoreline (border between Central and Southern region). Coastal marine waters tend to have higher PROP concentrations than off-coast. A west-to-east increase of pollution emerges also for particle-bound PROP concentration ( figure 2(b)). Some remarkable differences with figure 2(a) exist, the most evident being the Levantine basin, where the modelled particle-bound PROP concentrations typically exceed 1000 ng PROP /g microplastics but can even exceed 10 000 ng PROP /g microplastics , sometimes revealing strong geographic signatures (e.g. the cyan path-like pattern from the Nile Delta, in the Southeastern region, to the southern Turkish coast, Northeastern region).
While the mapping in figure 2(b) has no direct reference to particle counts, in figure 3(a) we show the spatial distribution (averaged over 2015-2016) of particles whose PROP concentrations exceeded the MPC. These highly polluted particles are distributed throughout the Mediterranean Sea, with the only exception of an area located close to the Strait of Gibraltar. The highest concentration of these particles is found in the Southeastern region, especially close to the Nile delta, and in the Gulf of Gabès and the adjacent coastline, in the Southern region. Aggregations of particles with high PROP concentrations are also found in the northern parts of the Aegean and Adriatic seas. To complement from the deep sea perspective the information about particle-bound PROP distributed in sea surface waters, figure 3(b) shows the spatial distribution of particle-bound PROP mass removed along with particle sinking during the whole simulation. PROP removal from surface through this pathway accounts for about 2.5 kg year −1 only at the whole basin scale, with the largest quantities sinking close to the Nile delta, and in general along the Levantine coasts. At the opposite side of the spectrum, particle-bound PROP sinking attains its minimum in a wide area extending approximately from the Ionian Sea to the Gulf of Sirte (at the crossroads of the Central, Northeastern, and Southern regions), as well as  in other smaller areas in the west coast of Corsica and close to Gibraltar.
The coupled model also interestingly permits to keep track of the chemical exchanges between seawater and particles (figure 4). Net adsorption/desorption of the focal PROP onto/from particles, integrated over the entire simulation, are distributed in vast, spatially clustered areas in each regional sea within the Mediterranean basin (see again figure 1(b)). High values of net adsorption are typically recorded along the coasts, especially in the Western Mediterranean ( figure 4(a)). By contrast, there are large areas in the Adriatic and Eastern Mediterranean Sea where particles tend to release the adsorbed PROP, like in the south-central part of the Adriatic Sea, the Aegean Sea and the northern Levantine Sea (Northeastern region). There is a wide portion of the basin, from the Ionian Sea to the Gulf of Sirte, where water-particle PROP fluxes appear to be weaker, similarly to what occurs in the Alboran Sea, albeit on a much smaller area. Although particles from all source types were clean of PROP upon release, net adsorption/desorption patterns differ when partitioning particles by their source type. Coastal particles ( figure 4(b)) show high desorption in the Adriatic and Eastern regions, yet PROP adsorption prevails in the Western and Southern regions of the Mediterranean basin. Riverine particles (figure 4(c)) appear to be responsible for important PROP releases throughout the Mediterranean basin, contributing not only to highdesorption areas (as those in the Southern Adriatic-Ionian Sea, and in the Levantine basin) but also to areas with weaker overall desorption, like south of the Balearic islands (Southwestern region). Particles released by fishing vessels (figure 4(d)) seem to mostly adsorb the simulated PROP, a result that is not surprising since fishing-related particles enter (clean, like all others) the sea from offshore locations, which are generally less polluted than coastal waters (see figure 2(a)). For this reason, and because of slower kinetics due to a smaller gradient between the particle-bound and seawater PROP concentrations (see equation (2)), particles released in fishing grounds are more likely to be chemically downgradient (i.e. they tend to adsorb) during their transport at sea. However, some areas with significant values of specific desorption can be identified in figure 4(d), mostly in the same regions as in figures 4(b) and (c).
Lagrangian tracking also allows to apportion microplastic-mediated chemical pollution, resulting from the coupling with Eulerian advection-diffusion of our target PROP at sea, to specific source types and regions of origin of microplastic particles. As an example, the origin of desorbing particles is elaborated for two focal areas, one located in the Ionian Sea, at the border between the Adriatic, the Central, and Northeastern regions ((a)-(e) in figure 5) and another between Crete and Cyprus (Northeastern region, (f)-(l)). In the Ionian area (figure 5(a)), the largest number of PROP-releasing particles was found to originate from the Central region ( figure 5(b)), and from coasts as source type (figure 5(c)). However, when considering the average mass of PROP desorbed per unit of microplastics' mass, particles from regions non-contiguous to the focal area, such as the Western and the Southeastern, also figured as significant carriers of pollution ( figure 5(d)). Among source types, riverine particles were those with the highest desorbed PROP mass-to-particle ratio in this area ( figure 5(e)). In the second focal area (figure 5(f)), the majority of desorbing particles originated from the Southeastern region (see the outstanding input from rivers, mainly driven by the Nile river and its delta, figure 5(h) and figure S1). Interestingly, despite not discernible from the pie plot in figure 5(g), all the seven regions appear as particle origins: not only a small number of particles came from the Southern and Central regions, but some even fewer from the Southwestern, Central, and Adriatic sub-basins. Figure 5(i) shows that particles sourced in the Northeastern and Southeastern regions are those that carry the most PROP into the area under exam. Conversely, particle source type (i.e. coastal, riverine, or fishingrelated) seems not to be significant in determining the amount of PROP desorbed by particles ( figure 5(j)). Additional analyses for other relevant areas are reported in SM, section S5.

Discussion
We present a first modelling framework designed to account for the chemical exchanges between microplastic particles and PROPs (Guerrini et al 2021) in the marine environment by coupling the dynamics of (Lagrangian) particles and (Eulerian) PROPs within the sea surface layer. The novel axis that this problem required us to explore is that plastic particles both (a) act as tracers passively transported by currents, and (b) chemically interact with seawater. From a technical perspective, we turned passive Lagrangian particles into active agents able to change the spatial distribution of PROPs via two-way interactions occurring at the particle scale. Although we are dealing with abiotic materials, some links exist between our approach and Lagrangian biophysical models, where the perception of relevant environmental variables can directly influence the trajectory of particles (see e.g. Melià et al (2013), Lett et al (2020), about the modelling of fish larvae, or of fish otoliths, Zenimoto et al 2011, Calò et al 2018.
In our model, plastic particles were assumed to be unpolluted at release: consequently, they take up a target PROP from the surrounding (polluted) environment. Interestingly, maps of average particlebound PROP concentration (figure 2(b)) do not exactly match those of average distribution of such PROP in seawater ( figure 2(a)). Particle history (including source type and location, other than journey at sea) influences the particle-bound concentration of the target PROP, as this depends on the chemical gradients experienced by the particle during transport, as well as on the consequent PROP exchanges. When quantifying them, PROP adsorption on particles prevails throughout the Mediterranean basin ( figure 4(a)), as expected because of the initial conditions at release. Contrasting figure 4(a) to figure S2(c), we observe that the target PROP tends to partition on particles mostly (yet not exclusively) close to source locations and especially along the coasts, from where 80% of the simulated particles is released (coastal and riverine sources). After release, transport of particles by surface currents contributes to exporting the PROP from coastal areas to high seas, reinforcing findings by De Frond et al (2019). This is the case, for example, of the desorption area in the Northeastern region (figures 4(a) and 5(f)) identified by our simulations, which can be attributed to particles released from the densely populated areas along the Nile River ( figure 5(g)). This meridional pathway, also visible in figure 2(b), is an occasional feature of the surface circulation of the Mediterranean Sea that occurred in the simulated period (see figure  4 in Millot and Taupier-Letage 2005). Nonetheless, this example is notable not only because the number of particles and the amount of PROP released in the Southeastern region are the highest in the whole basin (figure S1), but also because our simulations (coherently with the-scarce-field data available, see e.g. Berrojalbiz et al 2011) identify such area as highly polluted (figure 2(a)). South-eastern Mediterranean coastal waters are also a recurrent gateway for those particles that exchanged the most PROP during the simulation (figure S6, SM).
Our approach reveals that large-scale transport of particles and bounded PROP can even occur between distant sub-basins, as also visible in the focal area examined in figure 5(a), where some highly desorbing particles originated from the Western region (figure 4(d)) where particles tend to accumulate PROP ( figure 4(a)). In particular, the likelihood of microplastics to serve as long-range vectors of PROPs seems quite dependent on their source type. We found that riverine particles are the most significant in conveying the target PROP in the Mediterranean basin, largely contributing to the vastest desorption areas (figure 5(c)); moreover, they are responsible for the release of the highest amounts of the target PROP per gram of microplastics in some desorption areas (figures 5(e) and (l)). Since all particles were released clean, our results can once again be entirely ascribed to adsorption of PROP occurring during transport at sea, and not to processes occurring before release (e.g. industrial effluents within the river catchment). Figure 5(b) also highlights how specific features of the Mediterranean surface circulation can have important impacts on the dispersal of particles: Adriatic particles here account for about 25%, while in an adjacent location (see figure S3) desorbing particles sourced from this region make up 75%. The presence, in the southern Adriatic sub-basin, of counterclockwise gyres (Millot and Taupier-Letage 2005) possibly hinders the outflow of particles from the Adriatic to the Ionian Sea through the Strait of Otranto. Furthermore, our model can scout locations of possible ecotoxicological concern. For instance, figure 3(a) reveals that the areas with higher presence of particles whose PROP concentration exceeds the MPC (Verbruggen and van Herwijnen 2011), i.e. the eastern, coastal part of the Southern region and the Levantine coasts, coincide with areas where adsorption is most intense (figure 4(a)). There, strong adsorption could be caused by high seawater PROP concentrations, especially along the Levantine coasts (figure 2(a)). However, not always highly PROPcontaminated particles are found in highly polluted seawater-as shown in figure 3(a) for the Gulf of Hammamet, on the African side of the Strait of Sicily, which is also characterized by generally higher particle counts than the neighboring areas (see figure  S2(c)). Particles reaching that area could (a) have taken up PROP elsewhere, and/or (b) have persisted in the Gulf for longer, thus building up high particlebound PROP concentrations over time. Hypothesis (i) is enforced by previous results in Guerrini et al (2021), figure 4(c) and by the general circulation of the Mediterranean: particles entering the Southern region are likely to be conveyed by the Atlantic water stream from the Western regions (Pinardi and Masetti 2000). Surface circulation may also lend support to hypothesis (ii), though. In fact, currents passing the Strait of Sicily can also give origin to eddy-like circulation patterns south of the Strait of Sicily, as also shown by Pinardi and Masetti, potentially recirculating particles there and towards the coastline south of it. Figure 3(b) could hint instead to benthic areas where particle-mediated PROP contamination would occur. Interestingly, the PROP sunk with particles during our simulations is quite low, amounting to approximately 5 kg in 2015-2016. This could be an underestimated figure of the variety of PROPs that could be conveyed to depth by microplastics, as we considered phenanthrene only among the several types of pollutants that bind to plastics (Rochman 2015). Nonetheless, any conclusion regarding the benthic fate of pollution by PROPs has to be considered with caution, as we did not consider subsurface dynamics. In fact, the vertical trajectory of particles during sinking and settling could be far from a straight line from surface to seafloor, besides being potentially subject to ups and downs in the water column due to biofouling dynamics (e.g. Kooi et al 2017, Fischer et al 2021 or physical processes (de la Fuente et al 2021).
It is particularly worth reminding, and unfortunate, that little is known yet about the magnitude of the influxes of PROPs to the Mediterranean Sea: this knowledge gap heavily influences any quantitative assessment of modelled concentrations of PROPs, including ours, both in seawater and particle-bound. The major limitation faced by novel, holistic modelling approaches resides in lack of data. Our model would be potentially capable of more accurate projections if forced with actual influxes of PROPs, in terms of both quantities and source types. For instance, industrial activities (Tobiszewski and Namieśnik 2012) on land, as well as commercial shipping (UN/MAP 2017, UNCTAD 2020) and drilling rigs (Cordes et al 2016) offshore, should be encompassed as sources of PROPs. Furthermore, atmospheric deposition should also be accounted for as an input of marine pollutants to the Mediterranean Sea (Castro-Jiménez et al 2014). Finally, the interactions between microplastics and PROPs are strongly affected by polymer-chemical pairing (Tourinho et al 2019, Menéndez-Pedriza and Jaumot 2020); to exemplify the application of our approach, here we specifically focused on polyethylene particles and used phenanthrene as target PROP, also taking advantage of its high affinity and conservatively overestimating particle-bound PROP concentrations. Further modelling efforts could account for the variety of polymers and particle shapes among microplastics, influencing both their capacity of exchanging PROPs (Rochman 2015) and sinking behavior , and for the simultaneous presence of different PROPs at sea (e.g. Albaigés 2005).

Conclusions
Our modelling results help unravel the spatial patterns of particle-mediated adsorption and desorption of PROPs, and could thus represent a first quantitative evidence of the capacity of microplastics to act as vectors of pollutants, with focal application to the Mediterranean Sea. Microplastics confirm here their potential for being effective carriers of contaminants, able to exchange considerable amounts of PROPs with respect to their mass. Yet, particlemediated translocation of PROPs across the Mediterranean Sea (and to its depths) seems to be low at modelled (and realistic, see SM, section 2.3) concentrations of microplastics.
The Mediterranean Sea has microplastic concentrations among the highest in the world (Suaria et al 2016). To cut the effects of microplastics and marine pollution by PROPs in the long run, even a slight reduction in the amount of mismanaged plastic waste per capita produced in a country or higher efficiency in blocking/preventing microplastics from entering the sea would have a significant impact (see also Boucher and Billard 2020). This observation is made explicit in our proxy indicators for coastal and riverine microplastic input, where MPW production is a prominent factor in determining particle release (Tsiaras et al 2021). However, our results importantly suggest that curbing the input of microplastics by improving waste management might not be enough to tackle serious ecotoxicological risks, because even not abundant-yet highly pollutedmicroplastic particles can be found at sea as resulting from their actively exchanging transit through polluted areas.
To effectively fight the environmental threats posed by (micro)plastic pollution to the Mediterranean Sea and the global oceans, we thus ought not to forget about chemical pollution. In fact, due to high partitioning coefficients for several contaminantpolymer type combinations (O'Connor et al 2016), microplastics tend to concentrate contaminants by adsorbing them (León et al 2018), as also confirmed by our simulations. Future modelling efforts in this field should also explore pollutant fluxes to the seafloor and, in general, their distribution along the water column (as suggested by Jalón-Rojas et al 2019, Soto-Navarro et al 2020 and de la Fuente et al 2021 for microplastics), to extend risk assessment procedures from microplastics on the ocean surface (Everaert et al 2020) to PROPs and their related vertical distributions. Additionally, although the case study presented here is that of plastic-related chemical pollution in the Mediterranean Sea, our coupled approach is flexible enough to be applied to other marine regions, and to be complemented with other kinds of seawater-particles interactions (i.e. biofouling).
Healthier seas and oceans for future generations (United Nations 2020) need the implementation of comprehensive approaches, not only to identify and protect the mechanisms interconnecting marine ecosystems but also to understand the role played by the processes that are disrupting them. Quantitative methods, such as data-driven numerical modelling like ours, can shed some light on phenomena that are hard to investigate on the field. Approaches that couple physical risks for biota caused by microplastic pollution to chemical hazards related to PROPs, subject to potential biomagnification effects, are crucial to inform and prioritize actions on key pollution sources, a pivotal task for marine conservation.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).