Here we apply a model of ocean microplastic interaction with the base of the marine food web and biogeochemistry, embedded within an Earth System Model of intermediate complexity ( Figure 1 ). The microplastic model used in this study is modified from Kvale et al. (2020b) and fully described there. The marine biogeochemistry is described in Keller et al. (2012); Kvale et al. (2015b). The Earth System Model, the University of Victoria Earth System Climate Model (UVic ESCM) version 2.9 is described in Weaver et al. (2001); Meissner et al. (2003); Eby et al. (2009). Briefly, the UVic ESCM is a coarse-resolution (1.8° × 3.6° × 19 ocean depth layers) ocean-atmosphere-biosphere-cryosphere-geosphere model. The ocean levels increase in thickness from 50 m in the surface layer. Atmospheric processes are simplistically represented as a two dimensional energy moisture-balance model. The biogeochemistry resolves dissolved inorganic carbon and alkalinity, as well as oxygen, nitrate and phosphate ocean tracers. A seasonally cycling iron mask parametrizes iron limitation for primary producers (Keller et al., 2012). Phytoplankton are represented by three functional types: ordinary phytoplankton, calcifying phytoplankton, and nitrogen fixing diazotrophs (Keller et al., 2012; Kvale et al., 2015b). Calcifying phytoplankton and a single zooplankton class produce particulate inorganic carbon (carbonate) at a fixed ratio (4%) from particulate organic carbon production. This ratio was determined by model tuning in Kvale et al. (2015b) to achieve optimal fit against observed alkalinity and dissolved inorganic carbon (DIC) distributions (Key et al., 2004). Carbonate ballasts associated organic carbon and dissolves quasi-thermodynamically following the parametrisation of Aumont et al. (2003) but allowing for supersaturated dissolution (Milliman et al., 1999). The sediment model (Archer, 1996) is not activated for these simulations because water-sediment carbon interaction operates over timescales longer than those examined here.

The microplastics model represents three plastics partitions in the ocean using concentration units of particles per cubic metre (Kvale et al., 2020b). For simplicity, this model does not resolve particle size, type or breakdown. Unattached microplastic (MP) enters along coastlines and major shipping routes using a scaling against regional CO₂ pemissions to approximate population density and degree of industrialization. These unattached microplastic particles trace passively, but a fixed portion of the particles in each grid cell are prescribed a rise rate to approximate a quasi-buoyant body of microplastic in the ocean. Because size and type is not considered, all unattached microplastic is uniformly appealing and interacts equivalently with zooplankton and marine snow. Zooplankton ingest and then excrete the microplastic, which is then bound in a faecal pellet partition (MP_Z). Ingestion is calculated with a Holling II function which considers the grazing preference for microplastic over other types of food, and the total amount of food available (Kvale et al., 2020b). Faecal pellet microplastic sinks and is released back to the unattached plastic pool at the rate of detrital remineralisation. If the faecal pellet reaches the seafloor before returning to the unattached pool, a prescribed portion is considered lost from the ocean and the remainder returned to the unattached pool. The third partition is formed by the aggregation of the unattached partition with marine snow (MP _A ). Marine snowaggregates are produced as a fixed fraction from phytoplankton detritus formation, and microplastic aggregates with it using a Monod function to approximate a saturation threshold for the organic aggregate: microplastic ratio (Kvale et al., 2020b). Once bound in marine snow, this microplastic partition sinks at the rate of organic detritus and is released back to the unattached partition at the rate of detrital remineralisation. As with faecal pellet plastic, any marine snow bound plastic which reaches the seafloor either returns to the water column as unattached microplastic, or is considered lost to the system, with the return portion being prescribed.

Explicit three-dimensional modelling of biofouling vertical trajectories of several common microplastic sizes and types by others suggests that on average, biofouling has the potential to produce particle maximum concentrations within the upper 50 meters (Fischer et al., 2022), which would be completely within our model’s surface depth layer. The Kvale et al. (2020b) model represents this biofouling and defouling implicitly by biological uptake dynamics that include large losses to the temperature-dependent microbial fast recycling loop (Kvale et al., 2015a; Kvale et al., 2015b).

The microplastic model contains seven unconstrained parameters: the input rate of microplastic to the ocean, the fraction of unattached microplastic which rises, the fraction of phytoplankton detritus which aggregates into marine snow flocculates, the microplastic uptake coefficient for marine snow, the fractional return of microplastic at the seafloor, the zooplankton grazing preference for microplastic, and the microplastic substitution ratio for food in a zooplankton’s diet. These seven parameters were subjected to extensive sensitivity testing in Kvale et al. (2020b) and the reader is referred to the Supplemental Information accompanying that manuscript. While 14 of the parameter tests were identified as the most probable solution space, it was not possible to further narrow the candidate parameter sets. Thus, in that and subsequent publications (e.g., Kvale et al., 2020b; Kvale et al., 2021; Kvale and Oschlies, 2022) a set of three representative parameter combinations were used to describe robust features of the 14 member “most likely” ensemble.

In the Kvale et al. (2020b) model, the sinking rate of faecal pellets was not altered by the presence of microplastics. However, studies have shown that the sinking rate of faecal pellets can decrease when low density microplastics are included in their matrix (Cole et al., 2016; Coppock et al., 2019). To investigate the biogeochemical and ecological effects of slowed faecal pellet sinking in the presence of microplastic contamination, we use the experimentally derived linear relationship between increasing faecal pellet density (decreasing microplastic content) and increasing sinking rate (Cole et al., 2016):

The theoretical faecal pellet sinking rate (R _sink ) in meters per day depends upon a scaling (m; 411.4 meters cm³ day^-1 g^-1) multiplied by the faecal pellet density (ρ_pellet ) in g per cm^-3 and added to c (-432.0 meters per day), the y-intercept of Cole et al. (2016), their Figure 3A). This linear relationship was obtained by feeding polystyrene microbeads to C. helgolandicus, a species of copepod (Cole et al., 2016), and applies to purely organic as well as plastic-laden pellets.

We calculate the faecal pellet density for Equation 1 in each model grid cell from the density of organic material held in faecal pellets (ρ_{pellet 0}) and the density of microplastic (ρ_mp ) in g per cm^-3, each multiplied respectively by the proportion of mass of organic material (D) and the proportion of mass of microplastic in faecal pellets (MP_Z) for the grid cell:

Microplastic is simplistically assumed to have a density of 1.05 g cm^-3 as this is the density of polystyrene beads used by Cole et al. (2016). However, the base biogeochemical model does not explicitly represent organic particle density. It instead applies a linear rate of increase to all particle sinking rates (w_p (z)) following Berelson (2001). All organic detritus, including faecal pellets, sink at a rate from 14 meters per day at the ocean surface (w_{p 0} in Equation 3), increasing at a rate of 2 cm per day for each meter increase in model ocean depth (z) (to a maximum of 130 m d^-1 in the deepest parts of the ocean):

The motivation of this parametrisation is the observation of increased sinking rates with depth due to the loss of organic material by remineralisation as the particles settle (Berelson, 2001). When the organic material remineralizes, this leaves behind the denser siliceous or carbonate components of the pellet, thereby increasing the settling rates (Berelson, 2001). Thus in the base biogeochemical model, implicit particle densities for faecal pellets are assumed to increase from the surface to the deep ocean. This simplistic parametrisation is computationally efficient and produces a good agreement with observed nutrient, carbon, and oxygen distributions (Kvale et al., 2015b). Inserting the model pellet sinking rate range of 14 to 130 m d^-1 in Equation 1, we obtain the implicit model pellet density range of 1.084 (surface) - 1.366 (bottom) in g per cm^-3. Zooplankton faecal pellet densities vary widely in nature (e.g., 1.07 - 1.37 g cm^-3; Butler and Dam, 1994 and 1.7 g cm^-3; Ploug et al., 2008). We apply our calculated range of implicit faecal pellet densities to Equation 2 in our modelling framework, so that at each model depth level, and for each grid cell’s faecal pellet and faecal pellet microplastic content, a unique density ρ_pellet is obtained. This ρ_pellet is applied in Equation 1 to obtain the pellet sinking speed for each grid cell in the model.

Unit conversions are required in order to calculate the mass ratios of faecal pellets to microplastic in faecal pellets for each grid cell in Equation 2. Model units for microplastic are particles m^-3. Microplastic particle concentrations are converted to grams m^-3 using a conversion of 236000 tonnes = 51.2 ×10¹² particles (van Sebille et al., 2015), or 4.6 ×10^-3 grams per particle, as is already applied in other parts of the model (Kvale et al., 2020b). This conversion is similar to the 4.8 ×10^-3 grams per particle of the larger size (20.6 μm) polystyrene beads used by Cole et al. (2016). Mass conversion for the smaller sized (7.3 μm) polystyrene beads used by Cole et al. (2016) is 2 × 10^-4 grams per particle and results in overall denser faecal pellets for the same microplastic particle concentrations. Model units of organic faecal pellet material are mmol N m^-3. Nitrogen units are converted to grams of nitrogen using the molar mass 14 g N per mol. We use the average faecal pellet content (30 μg N: 1 mg faecal pellet) from Morales (1987) to convert to a mass in grams of pellets per grid cell.

To illustrate the effect of Eqns. 1 and 2 on Equation 3, we apply them to a published high surface plastic concentration simulation from Kvale et al. (2021), using model output of globally averaged organic detritus and microplastic concentrations with depth. This model simulation has been forced with rapidly increasing atmospheric CO₂ as well as rapidly increasing microplastic pollution rates, and therefore represents an extreme pollution and climate change scenario projected to the year 2100. When no microplastic is present, the calculated depth-dependent faecal pellet sinking rate is unchanged from the default model configuration value ( Figure 2 , black line). At increasing concentrations of microplastic, deviation of R _sink from R _{sink 0} becomes increasingly significant (see the additional slowing between years 2000 and 2100, dashed versus solid lines in Figure 2 ). The mass conversion of simulated microplastics to small or larger-size polystyrene can also exert a significant influence on calculated faecal pellet sinking rates (see blue versus red lines in Figure 2 ).

Note that the faecal pellet density calculation (Equation 2) occurs for each grid cell and thus the applied faecal pellet sinking rate depends strongly on the ratio of faecal pellet microplastic to faecal pellets in each grid cell. This parametrisation permits feedbacks, i.e. biogeochemical side effects of slowed surface pellet export due to microplastic contamination will affect both the microplastic and faecal pellet concentrations throughout the water column and therefore adjust sinking rates. This possibility of feedbacks is not represented in our offline calculations in Figure 2 , but is represented in our fully-coupled experimental framework which is described below.

We adapt the experimental framework previously used by Kvale et al. (2020b); Kvale et al. (2021); Kvale and Oschlies (2022). The UVic ESCM is first equilibrated at year 1765 boundary conditions, including agricultural greenhouse forcing and land ice. Historical CO₂ forcing is next applied from year 1765 to 1950. From 1950 to 2100 the model is then forced with both historical CO₂ forcing (to 2005) and a high-end CO₂ concentration projection (RCP8.5; Riahi et al., 2007; Meinshausen et al., 2011) ( Figure 3 ). Microplastic is introduced to the ocean model from year 1950 and the pollution rate increases 8.4% per year (Geyer et al., 2017). The total amount of microplastic added to the ocean is fractionally reduced from a calculated maximum potential pollution rate, which was tuned to fall close to the 275 million metric tonnes of plastic waste generated by coastal countries at year 2010 (Jambeck et al., 2015). We repeat this simulation sequence for three variants of the model; the first being the default, linearly increasing faecal pellet sinking rate version (CNTL) and the second (SLOW _big ) and third (SLOW _small ) being the density-modified slowed faecal pellet sinking rate version described by Equation3. The CNTL model includes zooplankton ingesting microplastic and therefore also simulates the corresponding biogeochemical response of ecosystems as described by Kvale et al. (2021). For SLOW _big we apply the mass conversion of large-size polystyrene microbeads (4.6 × 10^-3 grams per particle), and for SLOW _small we apply the mass conversion of small-size polystyrene microbeads (2.× 10^-4 grams per particle) from Cole et al. (2016). By comparing the model outputs of CNTL and SLOW, we can assess the biogeochemical impact of microplastic-slowed faecal pellet sinking rates in a warming climate. By comparing the model outputs of SLOW _big and SLOW _small we can assess the sensitivity of our model results to the assumed microplastic particle size. In a final step, we simulate a climate-change only scenario that does not include any microplastic pollution effects (CNTL _clim ). We do this in order to give climate warming context to the changes in ecosystems induced by microplastics pollution.

Microplastic model parameters remain poorly constrained. Therefore, in our exploration of the hypothesis that microplastic-slowed faecal pellet sinking impacts global ocean carbon cycling we apply our three model variants to 14 microplastic model parameter combinations. These combinations were obtained in Kvale et al. (2020b) using Latin Hypercube methodology (McKay et al., 1979). Three hundred unique microplastic parameter combinations, representing the whole of the possible parameter value combination space, were used to force the model between the years 1950-2100. About half of the simulations were found to be numerically unstable, and the survivors were further eliminated if the pollution inventory exceeded 9 × 10¹⁴ particles at the year 2010 (Kvale et al., 2020b). This cut-off was applied based on Eriksen et al. (2014) estimating 3 × 10¹² large-size microplastic particles in the upper ocean in the year 2010, plus a 1-2 orders of magnitude assumed sub-surface inventory. However, this estimate might be under-estimated by two or three orders of magnitude, as more recent surveys have found very high subsurface concentrations of small-size microplastics around gyres (Zhao et al., 2022). The model microplastic parameter combinations which remained numerically stable throughout the simulations, and which produced plausible global microplastic inventories, represent a range of microplastic model parameter values and vary with respect to their pollution rates, the fraction of unattached microplastic that rises, the fraction of phytoplankton detritus which aggregates into marine snow, the marine snow-microplastic uptake coefficient, the fraction of microplastic which is lost to the seafloor, the zooplankton grazing preference for microplastic, and the mass of microplastic which is assumed to substitute for the same mass of food in a zooplankton’s diet. A second Latin Hypercube sampling that includes updated elimination criteria is planned for the future, but is outside the scope of the present sensitivity study.

In the following section, we present the most robust biogeochemical responses across the model ensemble to microplastic slowing of faecal pellet sinking rates. We provide ensemble mean responses from the members of SLOW _small , SLOW _big , and CNTL, respectively, in order to provide a full-ensemble comparison. We also distinguish when model responses (SLOW - CNTL) do not agree in sign within the ensemble (indicating strong response dependency on parameter values), and calculate the standard deviation of model response (SLOW - CNTL) in each ensemble. If the standard deviation of the model response is smaller than the ensemble mean response (SLOW - CNTL) then we define the response as robust (indicating low response dependency on parameter values). By examining a high-pollution, high-CO₂ emissions scenario, we quantify the maximum multi-decadal perturbation potential of ocean biogeochemistry to the threat of microplastic pollution slowing faecal pellet sinking whilst undergoing strong climate change.

The ecosystem response to the slowing of faecal pellet sinking by microplastic determines the response in air-sea gas exchange of carbon. In our simulations, calcifying phytoplankton in the tropics benefit from surface nutrient additions caused by shoaled pellet remineralisation, and calcium carbonate (CaCO₃) production is generally enhanced ( Figure 10 ).

In the previous sections we demonstrate that the alteration of ecosystems, particle fluxes, and biogeochemical cycling that results from microplastic slowing of faecal pellet sinking has only a small impact on global air-sea carbon exchange. How significant are these changes regionally?

Similar to the previously described deoxygenation from microplastic ingestion by zooplankton (Kvale et al., 2021), we find that ecosystems display divergent responses, depending on whether they are limited by macronutrients or not. The schematic in Figure 15 summarises the ecosystem and carbon export effects that arise from microplastic slowing faecal pellet sinking. Slowed pellet sinking shoals remineralisation and releases both microplastics and nutrients closer to the ocean surface. In regions that are not limited by nutrients, such as the equatorial upwelling zones and high latitudes, additions of nutrients to the surface ocean do not affect primary producers. However, net detrital flux is increased by the erosion of top-down control on primary producers due to a positive feedback loop between zooplankton exposure to microplastic enhancing its ingestion, and slowing pellet sinking rates increasing their exposure. This feedback is prevented in upwelling systems where low-oxygen subsurface waters permit denitrification that limits nitrate availability downstream. In nutrient-limited environments, calcifying phytoplankton benefit from regional additions of nutrients that advect out of nutrient-replete regions as well as local enhancement of surface nutrients due to shoaled pellet remineralisation. Enhanced calcification reduces surface alkalinity and weakens the ocean carbon buffer, which leads to a minor global decrease in atmosphere-to-ocean carbon flux. Shoaled pellet remineralisation in the eastern tropical Pacific and Atlantic upwelling zones simultaneously increases the release of ocean carbon, causing a net shift to carbon held in the atmosphere. However, microplastic particle mass is a key assumption and larger microplastics can alter the balance between mechanisms, preventing the large-scale enhancement of calcification and actually increasing ocean carbon uptake. Thus, understanding ecosystem impacts of microplastics contamination in the real ocean will require careful characterization of both ecosystem structure as well as the microplastics themselves.

These simulated responses tentatively confirm the hypothesis derived from the Cole et al. (2016) study that a shoaling of pellet remineralisation can weaken the ocean uptake of atmospheric carbon, for smaller-sized microplastic particles. However, it is the adjustment of the whole ecosystem to the slowed pellet sinking that determines the response of air-sea gas exchange over the global ocean. Only the eastern tropical Pacific and Atlantic upwelling zones demonstrate a response similar to what is predicted, in that a shoaling of remineralisation enhances surface DIC and increases ocean carbon release. Elsewhere it is the adjustment of ecological competition between phytoplankton functional types and the reduction of top-down control that determines whether the ocean carbon sink will respond. This result serves as a reminder that extrapolation of experimental results from individual- to global-scale is not always appropriate, as has been noted in other biogeochemical modelling exercises, e.g. Taucher et al. (2022).

We report ecosystem changes due to microplastic-slowed pellet sinking that show both large and small variability across our model ensembles. Impacts in the Southern Ocean and the central tropical and sub-tropical ocean basins tend to be uniform in magnitude and sign across our SLOW _small ensemble, suggesting that these impacts (regional enhancement of organic particle fluxes, surface nutrients, NPP, and calcification) are robust and relatively quantifiable features of a microplastic-slowed pellet sinking regime. In contrast, corresponding impacts in the high latitude North Atlantic and North Pacific are more variable across the model ensemble and depend more on model parameter values. This robustness discrepancy might be explained bythe dominant mechanism being bottom-up (alleviation of nutrient limitation) in the Southern Ocean and subtropics but top-down (grazing suppression) in the high-latitude North Pacific and Atlantic. What this says about both SLOW model ensembles demonstrating poor agreement with respect to surface microplastics concentration trends is uncertain. Numerical instability of the SLOW _big ensemble was unfortunate and limits our discussion with respect to the impact of plastic mass assumptions, although we note the top-down impacts are apparently more robust than bottom-up impacts in this ensemble. Teasing apart the deeper mechanisms for parameter (in)sensitivities in regional responses is something we leave to the future.

Fortunately, our results suggest that over a century timescale the biological carbon pump responds by no more than about 1% in basin average to a microplastic-slowing of faecal pellet sinking rates. A weakening of remineralized carbon storage in the deep ocean is concerning over multi-century timescales and longer, and it is not clear from our simulations if sustained microplastic pollution would continue to weaken the biological pump over these longer timescales. Answering this question would require the application of microplastic particle breakdown, which our model does not resolve. Over the nearer term, the net air-sea gas exchange is hardly impacted, even when using the larger-sized microplastic particle assumption. This is an anticipated result because neither adjustment of particle fluxes nor the biological pump itself have been demonstrated to influence ocean carbon uptake over human timescales (Koeve et al., 2020; Wilson et al., 2022). Regionally, however, changes to ecosystem function can be significant and of the same order as climate warming impacts, with alteration of phytoplankton detrital fluxes exceeding 50% in the North Pacific and disproportionate impacts on low-productivity gyres. This finding is worrisome, given the amount of carbon released to the surface ocean as microplastic in our simulations is approximately 0.0006% of the amount released as atmospheric emissions over the same time period.

Beyond slowed pellet sinking, microplastics have been shown to have other adverse impacts on ocean ecosystems (recently reviewed by Galgani and Loiselle, 2021; Kvale, 2022). Microplastics can reduce zooplankton population size (Shore et al., 2021) by stunting egg production (Heindler et al., 2017) juvenile development (Jeong et al., 2017), and by increasing chance of death (Zhang et al., 2019). Other studies suggest microplastics can reduce growth rates of phytoplankton (Besseling et al., 2014; Sjollema et al., 2016). Microplastic aggregating with marine snow has been shown to increase the rate at which it sinks (Long et al., 2015; Porter et al., 2018). Additionally, microplastics also increase the rate at which marine snow aggregates (Michels et al., 2018). None of these impacts are considered in our modelling framework, but all potentially influence global biogeochemistry. Likewise, our simplistic approach neglects the full characterisation of microplastic polymer types found in the real ocean. Higher-density polymers such as PET have been found deeper in water column sampling (Galgani et al., 2022) and may therefore exert different ecological pressures to the low-density polystyrene we characterise our microplastics as. Polyethylene and polypropylene are less dense and more abundant than polystyrene (Pabortsava and Lampitt, 2020), with possibly a stronger impact on faecal pellet sinking. Likewise, because we do not account for microplastic particle size, our microplastic inventories and the relative exposure of the simulated food web may be considerably different to real-world exposures (Kooi and Koelmans, 2019; Kooi et al., 2021). Better quantification of very small microplastics (Zhao et al., 2022) to nanoplastic distribution and abundances can help to improve this aspect of our model. Neglecting size fractions also means we cannot account for variable egestion or zooplankton-mediated particle breakdown, which would impact the pellet density. Lastly, our model reduces the complex food webs of the real ocean to only three phytoplankton functional types and a single zooplankton grazer, and does not account for the buffering potential of flexible stoichiometry (Tanioka and Matsumoto, 2017) and community adaptation (Martiny et al., 2022), which could reduce the magnitude of our reported ecological impacts.

Despite the limitations of our modelling approach we conclude that marine microplastic pollution can influence ecosystems via regionally disparate mechanisms with globally coherent impacts on carbon cycling. These impacts are unlikely to be significant for carbon cycling over the next century, but warrant further investigation over longer timescales and with a larger variety of modelling approaches. Of particular interest to carbon cycling is the impact of environmentally-relevant microplastic contamination on biological calcification. More experimental data on this would be welcome. Other ecological impacts are potentially more significant and of the same order as climate change impacts. In terms of ecological impacts, quantifying the relative strength of multiple disruptive microplastic effects (zooplankton ingestion, slowed pellet sinking, and still others we do not consider here) will be crucial in quantifying net impact. Empirical demonstration that impacts on individuals scale to ecosystem level is expected to remain a challenge until in-situ data collection is improved. Established frameworks in the climate-carbon space, including metrics for assessment of non-CO₂ greenhouse gas climate impacts in terms of their CO₂ equivalence, can be usefully re-purposed for examining global-scale ecological impacts from plastics contamination. Such metrics that calculate CO₂ equivalence (or in our case, carbon equivalence of a potential impact) over a variety of timescales will be vital for defining the planetary boundary threat of plastics (Steffen et al., 2015) and will permit plastics researchers to exploit the rich literature on intolerable climate impacts (e.g. Petschel-Held et al., 1999; Kvale et al., 2012; Steffen et al., 2018) in defining such a boundary. Lastly, our finding of poor correlation between surface microplastic concentrations and ecological perturbation, and a high model dependence on microplastic model parameter values furthermore suggests that any confidence to be gained in estimates of biogeochemical impact must be solidly based in experimental work focussed on parameter estimation. We look forward to future developments in this area.