Research ArticleENVIRONMENTAL STUDIES

More than 1000 rivers account for 80% of global riverine plastic emissions into the ocean

See allHide authors and affiliations

Science Advances  30 Apr 2021:
Vol. 7, no. 18, eaaz5803
DOI: 10.1126/sciadv.aaz5803

Abstract

Plastic waste increasingly accumulates in the marine environment, but data on the distribution and quantification of riverine sources required for development of effective mitigation are limited. Our model approach includes geographically distributed data on plastic waste, land use, wind, precipitation, and rivers and calculates the probability for plastic waste to reach a river and subsequently the ocean. This probabilistic approach highlights regions that are likely to emit plastic into the ocean. We calibrated our model using recent field observations and show that emissions are distributed over more rivers than previously thought by up to two orders of magnitude. We estimate that more than 1000 rivers account for 80% of global annual emissions, which range between 0.8 million and 2.7 million metric tons per year, with small urban rivers among the most polluting. These high-resolution data allow for the focused development of mitigation strategies and technologies to reduce riverine plastic emissions.

INTRODUCTION

Plastic pollution in oceans and rivers is an emerging environmental hazard (1), and accumulation on riverbanks, deltas, coastlines (2), and the ocean surface (3) is rapidly increasing. Of all the plastics ever made to date, it was estimated that 60% has been discarded in landfills or in the natural environment (4). Plastic pollution poses threats on aquatic life, ecosystems, and human health (5, 6). Plastic litter also causes severe economic losses through damage to vessels and fishing gear, negative effects on the tourism industry, and increased shoreline cleaning efforts, adding up to US$1.26 billion per year for the Asian-Pacific Rim alone (7). Work on the origin and fate of plastic pollution in aquatic environments suggests that land-based plastics are one of the main sources of marine plastic pollution (8), either by direct emission from coastal zones (9) or by transport through rivers (10, 11). Riverine plastic transport remains understudied, especially in areas that are expected to contribute most to global plastic emissions into the ocean (12). A better understanding of pathways and transport mechanisms of plastic waste to and within rivers and the global distribution of riverine plastic emissions into the ocean is a prerequisite to developing effective prevention and collection strategies.

Previous attempts to estimate the distribution of global riverine emissions of plastic into the ocean (10, 11) relied on empirical indicators representative of waste generation inside a river basin. These assessments demonstrated a significant correlation between (micro)plastic concentration data collected by surface trawls in rivers, national statistics on mismanaged plastic waste (MPW) generation, and population density. For both studies, an empirical formulation was presented on the basis of this correlation, which was extrapolated to other rivers where data were not available. This resulted in predicted plastic (micro- and macroplastics combined) emissions of 1.15 million to 2.41 million metric tons (MT) per year (10) and 0.41 million to 4 million MT year–1 (11). These studies did not account for spatial distribution of plastic waste generation or climatological or geographical differences within river basins. According to these studies, the 10 largest emitting rivers contribute 50 to 61% and 88 to 94% to the total river emissions. Both models agreed on a disproportional contribution of Asian rivers to global plastic emissions. While these modeling efforts have provided a first approximation of the magnitude and spatial distribution of global riverine plastic emissions, they emphasized the scarcity of data on macroplastic contamination in freshwater ecosystems. Available measurements used for calibration of emission predictions were not always collected directly at the river mouths, and studies reported data on plastic contamination using varying units and methods, including surface trawling from boats or bridges (1315).

Sampling methods using surface net trawls for freshwater plastic contamination may be well suited for monitoring microplastic concentrations (size, <0.5 cm). However, insufficient sampled volumes limited by net opening width or pump outlet dimensions may result in the underestimation of macroplastics (several centimeters in size) (16) that account for most of the mass of plastic emissions (17). Instead, visual observations from bridges provide more consistent results for the quantification of floating macroplastic in rivers (18). In recent years, results from long-term visual counting campaigns for the quantification of floating macroplastic emissions from rivers of different continents have been made available (19). At a global scale, these studies provided observational evidence for the disproportional contribution of Asian rivers in plastic emissions predicted by numerical models (2024). Nevertheless, at a local scale, the studies reported discrepancies between observations and theoretical formulation (23), emphasizing the limitation of current models and the need for a revised formulation accounting for basin-scale geography, land use, and climate to more accurately estimate floating macroplastic emissions.

Here, we present a revised estimate of global riverine macroplastic emissions into the ocean using the most recent field observations on macroplastics and a newly developed, distributed probabilistic model to more accurately represent driving mechanisms of plastic transport (e.g., wind, runoff, and river discharge), differentiating between areas with different land use and terrain slope, and including plastic retention on land and within rivers. Microplastic transport is not included in this study. We derived probabilities for plastic waste to be transported from land to river and from river to sea from six different geographical indicators and generated a high-resolution (3 × 3–arc sec cells) global map of the probability for waste discarded on land to reach the ocean within a given year. This information combined with the most recent estimates of MPW generation on land (25) allowed us to estimate the annual emissions of plastic from rivers into the ocean. We calibrated and validated our model against 136 recent field observation data points (n = 52 for calibration and n = 84 for validation) of monthly riverine macroplastic transport from more than 67 rivers in 14 countries. We show how the consideration of transport probability for plastic within a river basin can highly increase or decrease the estimated emissions of the corresponding river into the ocean. At a global scale, this results in a considerably wider distribution of source points with large rivers contributing less to the total than expected, while urban rivers in South East Asia and West Africa are identified as the main hot spots for plastic emissions. We classified macroplastic-emitting rivers according to size, providing insight into which river class contains the highest number of rivers and the largest accumulative emissions. The classification and distribution of emission points provide a basis for development of mitigation strategies and technologies as well as a road map for upscaling existing mitigation technologies.

RESULTS

Study design

In this study, we calculate the probability for MPW generated inside a river basin to leak into aquatic environments. When combined with spatial data on MPW generation (25), our framework (Fig. 1) allows for the prediction of riverine plastic mass emissions ME into the ocean. Probabilities are derived from physical and environmental characteristics including precipitation, wind, terrain slope, land use, distance to river, river discharge, and distance to the ocean. Precipitation and wind are included as mobilizing forces and differentiate between climate zones, while terrain slope and land use are included to reflect the probability to reach a river and account for differences in landscapes from which MPW is generated. Distance to river and distance to ocean include the geographical location of MPW generation in relation to the nearest river and ocean. The probability of MPW to reach the nearest river depends on the landscape and the distance between the location of generation and the river. We conducted an expert elicitation to constrain model parameters and explored the impact of parameters with sensitivity analyses. We calibrated and validated our model against 136 field measurements of monthly emissions of floating macroplastics from 67 different rivers, collected between 2017 and 2020. A Monte Carlo and one-at-a-time (OAT) sensitivity analyses were performed, showing correlations of individual parameters with model output and field observations. On the basis of the ratio of residuals between 125 observed and modeled locations, a confidence interval was constructed. Model predicted emission points could range within a factor of 4 with a 68% confidence interval and a factor of 10 with a 95% confidence interval.

Fig. 1 Model framework.

Plastic emission in a river mouth ME is computed by accumulating of MPW multiplied with the probability of waste leaking into the ocean, P(E) within a river basin. P(E) is constructed with P(M), P(R), and P(O), which contain physical processes accountable for MPW transport.

Comparison with observations

A dataset of monthly averaged plastic transport near the river mouth was constructed from literature case studies and observational reports. This dataset was divided into a dataset used for calibration and one for validation (tables S1 and S2, respectively). Data that were collected and published before March 2019 were used for the calibration. Data published or made available to us after March 2019 were used for the validation. These studies use standardized methods to observe and quantify macroplastic transport according or comparable to published approaches (18, 26); see table S3 for details on observational and model-predicted riverine plastic emissions per month or per year.

Calibrated model results were compared with field observations, and a good order-of-magnitude relationship was demonstrated (coefficient of determination, r2 = 0.71, n = 51). All model predictions are within one order of magnitude from observations (the Pasig River is on the border of one order of magnitude) (Fig. 2 and fig. S1), except for the Kuantan River. The Kuantan River is considered an outlier, with observed concentrations an order of magnitude lower than estimated by the model; when the Kuantan River is included in the model, the coefficient of determination r2 is 0.61 (table S4).

Fig. 2 Observations compared with modeled data for floating macrolitter emissions per river.

Regression analysis carried out with 136 records from 67 different rivers of different sizes spread across the globe. The dataset was split into a calibration (n = 52) and a validation (n = 84) dataset. The coefficient of determination of the logarithmic regression, r2, is 0.71 for the calibration and 0.74 for the validation dataset. Symbols indicate midpoints of extrapolated measurements (MT month−1) on the x axis versus our best calibrated model prediction on the y axis. The horizontal whiskers indicate the upper and lower values reported for observational data (if published), and the vertical whiskers indicate the upper and lower value of the 68% confidence interval of model predictions. The dark blue symbols correspond to data points used for calibration, and light blue symbols represent the validation data points, while the symbol (triangle, circle, and square) indicates the continent from where the location originates. The logarithm of both the measurements and the model results is presented here. The dotted gray lines represent one-order-of-magnitude deviation from the x = y line in the middle. The Kuantan and Besos rivers (indicated in red) are outliers with more than one order of magnitude difference compared with observational results.

We validated (Fig. 2 and fig. S2) our model against 84 independent data points collected from literature on macroplastic observations or macroplastic correlated to microplastic observations. These data points originate from 51 rivers in six countries. We consider the Besos River (10 data points) an outlier because there are four weirs, which may act as a sink for plastics (27) directly upstream from the observation point (28). The remaining 74 validation data points are within one order of magnitude from the observed values, and our model predictions demonstrated a better correlation than during the calibration exercise (coefficient of determination, r2 = 0.74, n = 74). A separate validation graph (fig. S3) shows the results where 39 Japanese prefectures are merged into one single data point. This illustrates that the validation results are not biased by the inclusion of many points from one single country (r2 = 0.84).

Global distribution of riverine plastic emissions

Of the total 100,887 outlets of rivers and streams included in our model, we found that 31,904 locations emit plastic waste into the ocean, leaking in 1.0 (0.8 to 2.7) million MT into the marine environment in 2015. Rivers are included in the model if the annual average discharge is more than 0.1 m3 s−1 and counted as a plastic-emitting river if the annual plastic emissions are more than 0.1 MT year−1. Our model reveals that emissions are more widely distributed between contributing rivers with 1656 (range of 1348 to 1668 based on best calibrated lower and upper scenarios) rivers accountable for 80% of the global emissions against previously reported 47 and 5 rivers (Fig. 3A) (10, 11). In this study, we calculated a high-resolution distribution (3 × 3 arc sec) of probability P(E) for waste discarded on land to reach the ocean. P(E), with a global average of 0.4%, varied considerably between 0% for land-locked regions and up to 80% for coastal urban centers located near a river. When combined with distribution of waste generation on land, emission probabilities greatly increased the number of estimated riverine emission locations. This resulted in a considerably different ranking of the largest contributing rivers compared with previous assessments (top 50 rivers presented in table S5), as well as predicting the emergence of small rivers in the top ranking, for example, the Klang River in Malaysia.

Fig. 3 Global distribution of riverine plastic emission into the ocean.

(A) Contribution of plastic emission to the ocean (ME) (y axis) is plotted against the logarithm of the number of rivers accountable for that contribution (x axis), for previous studies and this study. (B) Distribution of 1656 rivers accountable for 80% of emissions over five discharge classes (x axis). Each river is represented by a dot. Within a discharge class, the position of a river (dot) is determined by the plastic emission (y axis). The boxes contain 50% (Q1 until Q3) of the data, and the solid horizontal line in the box is the median, while the dotted horizontal line represents the average emission per river within the discharge class.

On the basis of recent field observations and by considering probabilities of transport of plastic waste on land at high resolution within a river basin, we showed that land use, distance from waste generation to the nearest river, and coastline play a more important role than the size of the river basin itself. Hence, coastal cities associated with urban drainage and paved surfaces presented the highest emission probabilities, particularly in regions with high precipitation rates. On average, river basins with the dominant land use “artificial surfaces” are calculated to have a larger probability to emit plastic into the ocean than river basins with predominantly “cultivated land” (13 and 2%, respectively) and are observed and modeled to emit larger fractions of plastic waste into the ocean (15 and 3%, respectively) (see tables S6 and S7). To illustrate this, we compare the Ciliwung River, Indonesia and the Rhine River, Western Europe. The Ciliwung River basin on Java, which is much smaller than the Rhine River basin (respectively, 591 km2 versus 163,000 km2) and with less total generation of plastic waste (respectively, 19,590 MT year−1 versus 34,440 MT year−1), emits two orders of magnitude more plastic into the ocean (308 MT year−1 observed and 205 MT year−1 modeled for the Ciliwung River, and 3 MT year−1 observed and 5.4 MT year−1 modeled for the Rhine River). This difference is caused by the differences in geometry and climate between the two river basins and the spatial distribution of waste generation. In the Ciliwung River basin, waste is generated at 1 km from the river network on average and 29 km from the ocean. Waste generation in the Rhine River occurs, on average, at a much greater distance from the river network and the ocean with an average of 5 and 1021 km from the river network and the ocean, respectively. Moreover, the annual precipitation (29) in the Ciliwung River basin is more than 2.5 times larger than for the Rhine River basin (2445 mm year−1 against 950 mm year−1), further increasing the mobilization of plastic waste. The resulting average probability of waste to reach the river mouth for the Ciliwung River basin was 15.7% versus 0.04% for the Rhine.

We divided the 1656 rivers accountable for 80% of emissions over five river discharge, Q (m3 s−1), classes (Figs. 3B and 4A). The first class contains rivers with Q < 10 m3 s−1; the second class with 10 m3 s−1 < Q < 100 m3 s−1; the third class with 100 m3 s−1 < Q < 1000 m3 s−1; the fourth class with 1000 m3 s−1 < Q < 10,000 m3 s−1; and the fifth class with Q > 10,000 m3 s−1. We found that the 830 rivers in the first class combined account for 30% of global emissions. Middle-sized rivers, 582 and 211 in classes 2 and 3, respectively, combined account for 47%. Both in number (27 and 6 rivers in classes 4 and 5, respectively) and in combined emissions (2 and 1%, respectively), the large rivers account for a relatively small fraction. The remaining 20% of emissions is divided over 30,248 rivers of varying size and low (<92 MT year−1) emissions per river. Our results therefore suggest that focusing on implementing mitigation measures in small- to medium-sized rivers already could considerably reduce plastic emissions.

Fig. 4 Global emissions of plastic into the ocean.

(A) The geospatial distribution of plastic entering the ocean through rivers. The 1656 rivers accountable for 80% of the total influx are presented. The gray shading indicates the probability for plastic entering the ocean [P(E)] on a 10 × 10–km resolution. (B) Total emitted plastic into the ocean ME per country divided by the national generation of MPW, globally ranging between 0 and 20%. (C) Total emitted plastic into the ocean ME (MT year−1) per country.

Predicting national emissions and potential for plastic waste leakage into the ocean

We estimated that 1.5% (range, 1.2 to 4.0%) of the 67.5 million MT (25) of total globally generated MPW enters the ocean within a year. However, on a national level, the fraction of discarded waste entering the ocean differs considerably between countries (Fig. 4B). Our results indicate that countries with a relatively small land surface area compared to the length of their coastline and with high precipitation rates are more likely to emit ocean plastics (table S8). Particularly, for areas in the Caribbean such as the Dominican Republic and tropical archipelagos such as Indonesia or the Philippines, this results in a higher ratio of discarded plastic waste leaking into the ocean, respectively, 3.2, 6.8, and 8.8%. The plastic emissions of these countries are therefore disproportionally higher compared to countries with similar MPW concentrations but different geographical and climatological conditions. For example, Malaysia generates more than 10 times less MPW than China (0.8 million MT year−1 in Malaysia against 12.8 million MT year−1 in China); however, the fraction of total plastic waste reaching the ocean is 9.0% for Malaysia and only 0.6% for China. The largest contributing country estimated by our model was the Philippines with 4820 rivers emitting 356,371 MT year−1 (8.8% of the total generated MPW in the country), followed by India with 126,513 MT year−1 (1.0% of total generated MPW through 1169 rivers), Malaysia with 73,098 MT year−1 through 1070 rivers, and China with 70,707 MT year−1 through 1309 rivers (see Table 1 and Fig. 4C).

Table 1 Country statistics.

Top 20 countries ranked according to annual plastic emission ME into the ocean, as calculated in this study. The third column contains the annual MPW generated in each country. The fourth column contains the fraction (%) of MPW reaching the ocean (calculated by dividing national ME by MPW) within a year. The fifth column contains the country averaged probability for a plastic particle to reach the ocean within a year, P(E). This sixth column contains the number of rivers accountable for national emission ME, and the last column holds the number of rivers for a country that contribute to the global 80% riverine plastic emission (emitted by 1656 rivers in total).

View this table:

DISCUSSION

Our study predicts that riverine plastic emissions into the ocean are distributed across a much larger number of rivers than reported in previous modeling studies. The number of rivers responsible for 80% of global emissions (1656 in this study) is one to two orders of magnitude higher than previously reported [47 rivers (10) and 5 rivers (11)]. An important difference is that in previous studies, MPW was lumped within a river basin, leading to disproportionally high predictions of plastic emissions for large rivers while smaller rivers may have been underestimated. In this study, we considered the spatial variability of MPW generation within a river basin and introduced climate and terrain characteristics to differentiate the probability for waste to leak into rivers and subsequently the ocean. Therefore, MPW near a river and near the coast has a relatively high probability of entering the ocean, while MPW far upstream in a basin has a lower probability of entering the ocean. By considering these parameters, relatively small yet polluted river basins contribute proportionally more compared to large river basins with equal amounts of MPW generation within the river basins. Cities such as Jakarta and Manila are drained by relatively small rivers, yet observations and our model suggest that these rivers contribute more than rivers such as the Rhine or the Seine, for which the MPW generation is similar yet located further upstream.

The results from this study are important for the prioritization and implementation of mitigation strategies. The large number of emission points predicted by our model calls for a global approach to prevent, reduce, and collect macroplastic waste in aquatic environments instead of focusing on just several rivers. Furthermore, our results suggest that small- and medium-sized rivers account for a substantial fraction of global emissions. Predicted emissions presented in this study suggest that, besides the annual plastic emissions into the ocean, most plastic waste (98.5%) remains entrapped in terrestrial environments where it accumulates and progressively pollutes inland (aquatic) ecosystems. As most MPW is generated and remains on land, prevention and mitigation regulations for land-based waste reduction, collection, and processing as well as cleanups will naturally yield the largest impact on reducing the emissions of plastic into the ocean.

Understanding the total annual global riverine plastic emissions into the oceans is an important input for mass balance exercises and mapping the severity and fate of macroplastic pollution in the ocean. We calculated the annual global emissions to be between 0.8 million and 2.7 million MT. This in the same order of magnitude as previous river emission assessments, which estimated 1.15 million to 2.41 million MT (10) and 0.41 million to 4 million MT (11) for global riverine plastic emissions. However, a wider distribution of emission points in this study led to a new ranking of top contributing rivers, where the Pasig in the Philippines is now the most polluting river. It is good to note that, in comparison with the observational data from the Pasig river, the model overestimates the plastic emissions for this river; however, measurements were taken at the end of an unusually dry season (2019). This highlights the importance of including more detailed temporal variations in future model and observational studies, as well as the fact that model predictions could vary up to a factor of 10 (with 95% confidence) from field observations. Confidence intervals (95%) are therefore included in the list of the 50 most plastic-emitting rivers (table S5). The Yangtze River, which was previously estimated as the highest plastic-emitting river (10, 11), is now ranked 64th by our model. The Yangtze catchment is one of the biggest river basins with a very high total amount of MPW generation. However, the distance from MPW generation to the river and to the ocean is considerable as well. Therefore, according to our model, only a relatively small fraction of MPW reaches the Yangtze River and even a smaller fraction subsequently reaches the ocean. It is important to note that we calibrated our model against visual observations of macroplastics (>0.5 cm in size) and, therefore, we are not considering microplastic transport. Global riverine microplastic emissions are estimated to be several orders of magnitude lower (17) than our macroplastic emission estimate (47,000 MT/year, for a recent past scenario). Although plastic observations are extrapolated to the entire water column, our model does not include riverbed transport of plastic waste. Hence, our global riverine emission estimate can be considered conservative in the sense that it likely underestimates the total flux. We note that our estimated range for emissions in 2015 is one order of magnitude lower than a previous prediction for plastic waste inputs from land into the ocean (9) for 2010 (range, 4.8 million to 12.7 million MT year–1). This study did not specify a transport mechanism and includes all emissions into the ocean and not only riverine emissions, which emphasizes the uncertainty related to estimating plastic waste generation and emissions as well as the need for additional ground truth data.

Previous studies (10, 11) on global river emissions of plastic in the ocean were mainly calibrated against data collected in European and North American rivers. Following the recommendations from these studies, we included more data from South East Asian rivers to refine our model predictions. The difference between observed and calibrated modeled emissions is within one order of magnitude for 51 of 52 observational data points. Our model validation against 84 observed emissions yields 74 of 84 data points within one order of magnitude, with a clear explanation (four weirs upstream) for the 10 data points outside one order of magnitude. Given the uncertainty in observational accuracy and MPW data, we consider this an acceptable result and a major improvement compared to the performance of previous models. This study is limited to monthly average and annual emissions intended for the quantification of global riverine plastic transport and river-to-river comparisons. We expect temporal variations in discharge, and especially floods, to have a large impact on macroplastic mobilization and transport, as was found for microplastics (30). Therefore, future studies should include a higher resolution for temporal hydrological variations, aimed at better accounting for extreme events such as floods and quantifying their contribution to emissions. The model parameters chosen for this study are based on expert elicitation and calibration on field observations. More research and data are required to improve and validate the established relationships in this study. It is important to note that this study does not differentiate between types and characteristics of plastic waste. Mobilization, transportation likelihood, and buoyancy may be influenced by plastic particle properties such as shape, weight, and density. Therefore, the transport of plastic of different types and sizes should be differentiated in future assessments. Our global model does not include changes in national waste management policies nor the contributions of the informal recovery sector. It is good to note that the input MPW has a great impact on the results. Future research into plastic waste generation, practices, and distribution could substantially improve the accuracy of the model output. We also do not consider the presence of regulating structures in rivers such as dams, weirs, or trash racks and local extraction efforts, which have been reported to act as sinks for macroplastics (27). Depending on the operational characteristics, artificial structures in rivers, such as dams, may have a substantial impact on the downstream transport of plastics (10). We acknowledge the need for local modeling and observational studies to better address local conditions. The sensitivity analysis shows that the model is sensitive to changes in the probability for plastic waste to be transported 1 km downstream in a river. Future fieldwork aimed at measuring the trajectory (distance traveled downstream in time) of plastic particles in rivers is recommended. The uncertainty in parameter values should be minimized by conducting extensive monitoring campaigns on plastic mobilization and transport behavior rather than extensive calibration. Population densities, waste management practices, and consumption patterns are subject to change, leading to a varying generation of MPW (25). Ongoing efforts to improve global datasets on land cover, precipitation, and elevation continue to deliver more accurate input datasets. Our probabilistic modeling approach and framework allow for the inclusion of these improved datasets and benefit from parameterizations derived from local models with high-resolution temporal and spatial data on plastic transport and hydrology.

Our results include a global dataset of 31,904 locations representing river mouths and their estimated emissions. It is important to consider that emission values of individual rivers predicted by our best calibrated model may vary around a factor of 4 (68% confidence range) to a factor of 10 (95% confidence interval) from field observations. These data will be publicly available for researchers, policy-makers, and citizens to identify and address the nearest polluting river. Insights from this study should be combined with local field studies and local context from various scientific fields to effectively address plastic pollution in rivers and design and implement mitigation measures and strategies.

MATERIALS AND METHODS

Model formulation

The probability P(E) for plastic waste, discarded on land, to be emitted into the ocean is constructed from the probability of intersection of three events: M (mobilization on land), R (transport from land to a river), and O (transport from the river to the ocean)P(E)=P(M  R  O)=P(M)×P(R)×P(O)(1)

For each 3 × 3–arc sec grid cell, the amount of plastic waste leaking into the ocean is therefore calculated by multiplying the probability P(E) with the total amount of generated MPW mass (kg year−1) within the cell. The total annual emission ME of plastic into the ocean from a river is then computed by accumulating this product for all n grid cells contained in the river basinME=nMPW×P(E)(2)

Similarly to sediment (31) and debris (32), plastic waste may be mobilized during events of rainfall (33) where surface runoff accumulates and entrains plastic waste over the surface. Larger rainfall events may trigger buoyant transport of plastic litter. The amount of runoff, which is a fraction of rainfall, can vary depending on the land use, characteristics of the rain event, temporal variability of rainfall, and its distribution over the year. These variations can have a substantial impact on runoff and, therefore, on plastic transport triggered by rainfall. The relationship between runoff and plastic transportation, including hydrological thresholds, can be different depending on river basin characteristics. In this study, we focus on spatial variability and differentiation. For this study, we assume that annual wind and precipitation data are a good proxy for transport and mobilization probability. Global precipitation analysis suggests that areas with higher total annual precipitation also experience more frequent and intense rainfall events (34). In areas with little or no rainfall, wind can still mobilize and transport littered plastic waste on land, particularly from open-air landfills (35). Previous studies have reported positive correlations between both wind and rainfall with riverine plastic observations (28). If plastic mobilized by wind reaches a river, drain, or creek, then the flow dynamics will further transport plastic downstream. We do not include a directional component for wind in our approach. For this global study, we aim to differentiate between regions with little or high precipitation and wind and approximate the complex plastic mobilization mechanisms; we do this by assuming a linear relationship between annual total rainfall and average wind speed and their corresponding mobilization probability. These relationships are described in table S9. In our framework, we consider that plastic waste can be mobilized through both events of precipitation and wind. Hence, the probability of mobilization P(M) can be formulated from the union probability of precipitation event P and wind event WP(M)=P(P  W)=P(P)+P(W)(3)

Probabilities of mobilization by precipitation and wind are linearly ranging from 0% (respectively, no rain or no wind) to 100% corresponding to an upper threshold (see table S9). For probability of mobilization by wind, we consider the maximum monthly average wind speed (m s−1). The upper threshold for total mobilization was set at 32.7 m s−1, which equals to Beaufort 12 (i.e., under hurricane conditions, 100% of littered waste is mobilized). The upper threshold for probability of mobilization by rain was determined during the model calibration exercise presented later in Model calibration, considering the annual rainfall. Data for monthly averaged wind speed and annual rainfall were sourced from global 30–arc sec datasets distributed by WorldClim2 (29).

For the mobilized fraction of plastic waste, we compute the probability to reach the nearest river. The river network in our model contains the annual average discharge (m3 s−1) on a 3 × 3–arc sec spatial resolution and was derived by accumulating annual average 0.5° × 0.5° runoff between 2005 and 2014 (mm year−1) (36) by a nearly global flow direction grid (37). Cells with a discharge higher than 0.1 m3 s−1 are considered rivers (38). The shortest downslope distance DLand (km) from each grid cell to the nearest location in the river network is calculated on the basis of flow direction data. The result of this calculation step is a 3 × 3–arc sec grid containing distance to the nearest river.

The mobilization probability and the travel distance to the nearest river are now determined, yet the probability to reach a river also depends on the type of landscape between the location of mobilization and the nearest river. The objective here is to appreciate differences between, for example, flat and mountainous areas or rural and urban areas. Similarly to Chezy’s formula (39) and the rational method (40) in hydrology, we introduce a “roughness” coefficient based on land use classification and the terrain slope to approximate the complex mechanisms of overland transport of plastic waste. Here, roughness of a cell is defined as the transport probability to the nearest downslope cell. Rougher cells have a lower transport probability, while smoother cells have a higher transport probability. For example, plastic waste will more likely be transported by wind or rain on paved surface than in dense vegetation (33, 41). Higher-plastic observations are reported near the outlets of storm drainage, associated with urbanized land use, while areas further away from infrastructure have lower reported plastic litter concentrations (41). Vegetated areas, such as forests, provide more shelter and entanglement than bare areas, reducing the transport and mobilization probability by wind. Differences in roughness within a land use class including directional components may have a substantial impact on local scale but are not included in the global spatial resolution of the land use data used in this study. Furthermore, the average terrain slope (%) is known to increase erosion rates and sediment transport over land (42). Hence, the probability of transport to a river will naturally increase with terrain slope. The land use data were sourced from 30 × 30–arc sec classification distributed by GLC2000 (43) and contains 23 different land use classes (table S10), ranging from tree cover, shrub land, and cultivated areas to artificial and urban areas. The terrain slope was calculated from the 3 × 3–arc sec digital elevation model provided by HydroSHEDS (37).

We derive the roughness of each cell from land use and terrain slope and compute the average probability from the initial emitting grid cell to the nearest river cell, according to equations in table S9. As roughness is cumulated on the downslope path, the resulting probability to reach a river is exponentially decreasing with distance to river DLand. The probability of transport to a river is formulated as followsP(R)=(i=1nνi×(ε×si+τ)n)DLand(4)where νi is the probability for plastic to be transported to the next grid cell associated to land use (see classification in table S10) of grid cell i, si is the percent slope of cell i, ε is the coefficient that determines the increasing probability with increasing slope, τ is the minimum threshold probability for areas where the slope is zero (table S9), and n is the number of cells from origin to the nearest river cell.

By analogy to the transport of leaves (44) and wooden debris (45) by rivers, the probability in our model for plastic introduced in rivers, to reach the ocean, increases with river discharge and decreases with distance to ocean. Rivers with a higher Strahler (46) stream order (SO) have a larger cross section (47) and, therefore, on average, less friction (39), decreasing the likelihood for floating macroplastic to be intercepted. Therefore, for each river grid cell, we compute the distance DRiver to the ocean, the Strahler stream order, and the annual river discharge (m3 s−1). The probability for transport into the ocean is calculated as followsP(O)=(i=1n(θ×SOι+ι)×(κ×Qi+μ)n)DRiver(5)where θi is the probability related to Strahler stream order for cell i and ι is the lower-threshold rivers with the lowest stream order. Qi is the river discharge at cell i, κ and μ are the probability coefficient for increasing discharge and lower threshold probability, respectively (table S9), and n is the number of cells from the river entry point to the ocean. An example of the different steps leading to the calculation of probability of emission P(E) is provided in Fig. 5.

Fig. 5 Probability maps.

(A) The Meycuayan and Tullahan river basins and river network in Manila, Philippines. (B) The distance (km) from a 3 × 3–arc sec grid cell toward the nearest river. (C) The distance (km) from each grid cell to the ocean, through the river network. (D) The probability for a grid cell to emit plastic waste into the ocean P(E), Eq. 1, for a given year, ranging from 0 to 5% for areas further away from a river up to 0.8% for areas near a river and near the coast.

Observation-based estimates

We calibrated and validated the river plastic emission model using observation-based estimates. We used three types of data: (i) visual observations of floating macroplastics from own collection, (ii) published datasets of visual observations of floating macroplastics, and (iii) a published dataset of collected floating plastic debris. The visual observations were done according to the method developed for the Riverine and Marine floating macro litter Monitoring and Modeling of Environmental Loading. For this method, observations are conducted from bridges near river mouths. During each measurement, all floating plastic items are counted for a certain duration. The data are normalized over time by expressing the data in floating plastic items per hour (items hour–1). In general, the minimum observed item size is 0.5 cm, and items are visible in the upper 10 cm of the water column. We collected 136 data points from 67 rivers in 14 countries. Data that were collected and published before March 2019 were used for calibration, and data published after March 2019 were used to validate the model. Hence, we used 52 data points for calibration and 84 data points for validation of the model. The published field observations used for both calibration and validation were measured in rivers that have mutually different characteristics regarding total basin area, average land use, rainfall, and MPW generation (tables S6 and S7) and were collected on three different continents (fig. S4). Statistics on characteristics (minimum, median, maximum, first quartile, and third quartile) of the full observational dataset, the calibration dataset, and the validation dataset are listed in table S11. This table shows that the calibration and validation datasets have comparable statistics to the total dataset. We consider the data used for both calibration and validation representative as they are collected in different types of rivers and geographies while the total datasets have similar statistics. We use time of publication as criterion to include data in the calibration or validation. An overview of the observational data used in this study is presented in table S1 (calibration) and table S2 (validation) (4853).

We used a three-step approach to estimate the total plastic mass transport from the floating plastic observations. In the first step, plastic item flux p (items/hour) is converted from items/hour to mass flux using the mean mass per item c (kg per item). For each river, we either used the mean mass per item based on available data or used a lower and upper estimate based on the range available in the literature (0.002 to 0.019 kg per item). Second, we estimated the plastic mass flux in the upper 1.5 m of the water column using the floating fraction. For this fraction, we either used the value from the field observations or used a lower and upper estimate based on the global range available (0.6 to 0.95%). Last, we estimated the total mass transport over the entire water column using the fraction of transport in the upper 1.5 m. Here, a range was used with lower and upper values of 66 and 79% based on measurements in the Danube (20). Each step was performed using a range of values, yielding a lower and upper observation–based estimate. We used the middle value for model calibration and validation and use the range to quantify the uncertainty. Note that for the data already expressed as mass transport (e.g., Jones Falls and Japanese rivers), the first step was omitted. Data were collected using visual counting measurements of floating macroplastic litter from bridges (18, 26). This was converted into mass flux (MT year−1) using the following equationMobs=p×mp×c(6)

With observed floating plastic transport p (items hour−1), mp mean mass per plastic item (kg per item), and conversion factor c to account for plastics at deeper layers. We used both monthly and annual estimations in the comparison with the model results. Variables mp and c were measured at each river through net sampling at the same location as the visual counting measurements. If these measurements were not available, then we used the global or regional average values (table S1).

Expert elicitation and initial parameter bandwidths

To constrain our model parameters, an expert survey was conducted during the European Geosciences Union General Assembly, April 2019, in Vienna, with a panel of 24 geoscientists. In addition, the same survey was conducted online in 2020 to obtain a more global perspective from 36 geoscientists and other professionals active in plastic research. In total, this resulted in a panel of 60 experts from 39 organizations in 23 different countries across five continents (see fig. S5 for the global distribution of the expert panel). The advantage of benefitting from the intuitive experience of experts to assess complex modeling problems has been reported for hydrology (54) and ecology (55). Here, a series of seven questions related to the probability of plastic waste transport over land and through rivers were asked to individual experts. The questions are presented in table S12, while the individual responses are given in table S13. From this elicitation exercise, we calculated the average and SD of returned values for each question (table S14). These data determined the initial bandwidth for our parameters during the model calibration (i.e., while varying our model parameters when comparing with measurements, the resulting probability should remain in the range determined by experts elicited for this study, avoiding unreasonable parameter values).

The bandwidth for probabilistic parameters in our model can range between 0 and 100%. The parameter bandwidths presented in table S14 have been reduced on average by 46%. Sensitivity analyses, calibration, and validation further confine and support the parameter bandwidths and result in a parameter set that is in line with expert expectations and is optimized for the highest correlation (coefficient of determination) with field observations and the lowest ratio between observed and predicted plastic mass.

Model calibration

Our model predicts annual plastic emissions, which are scaled by monthly average discharge, to distribute these emissions over 12 months. As is common for complex environmental models, we calibrated our model, on the basis of output from the expert elicitation and sensitivity analysis, in an iterative way (56). First, we ran a version of the model to match the average values reported by the expert elicitation exercise. We evaluated the model performance by calculating the regression coefficient r2 between the logarithms of measured and modeled monthly averaged emissions and calculated the ratio between the sum of all observed and modeled plastic emissions. Under these conditions, the model-estimated emissions appeared higher than observations. We initially decreased the probability for plastic waste to be transported from land to a river cell P(R) by progressively increasing the roughness related to land use, as introduced in Eq. 4.

Second, the model overestimated emissions of rivers where precipitation was relatively higher than other rivers, when compared to observations. We adjusted our model results by decreasing the probability of mobilization P(M) induced by precipitation, as introduced in Eq. 3 based on results from the sensitivity analysis.

Third, the emissions of river basins where the generation of MPW is occurring further away from the mouth were underestimated (e.g., the Motagua in Guatemala and the Seine in France). Therefore, we altered our model predictions by increasing the probability of transport from river entry to ocean P(O), as presented in Eq. 5.

This model calibration exercise resulted in eight global iterations that are presented in table S4, showing the score model versus measurement per iteration, for the different parameters considered by our model. Our best calibrated scenario returned a regression coefficient of determination r2 = 0.71 between modeled and measured logarithm of monthly average emissions per river and with 51 data points modeled within one order of magnitude (fig. S1) from measurements. For this parameter set, we found that the sum of 51 predicted emissions exceeds the sum of observed emissions by only 6%, which we consider acceptable, in combination with the relatively high r2.

Workflow sensitivity analyses

Sensitivity analyses were performed to assess the sensitivity of the output of the model (tons of macroplastic emission per year per river mouth, ME) to parameter variations of all parameters in model Eqs. 4 and 5. Both a global sensitivity analysis and a local sensitivity analysis were conducted to assess the impact of interactions between varying parameters and the impact of varying an individual parameter. The global, all-at-a-time (AAT) sensitivity analysis provides insight into the interactions between input parameters and is implemented in this study using Monte Carlo techniques. The local, OAT sensitivity analysis provides insight into the response of the model to (extreme) variations of an individual parameter, while all other parameters are fixed. Parameters for the AAT method were randomly sampled from uniform distributions and cover a larger bandwidth than the boundaries provided by expert elicitation to explore a broader parameter space (table S15). For the AAT analysis, three river basins were selected to balance computational time and representability. For these basins, the model output of 100 parameter sets was obtained. The three river basins used for the simulation runs are the Tullahan, Motagua, and Jones Falls. These river basins have different annual rainfall, surface area, dominant land use, and average distances to the river (within the basin) and to the coast (within a river) (characteristics listed in tables S6 and S7). For the OAT analysis, each parameter was decreased and increased with an order of magnitude, while the remaining parameters remained constant (table S16), resulting in 16 simulation runs (2 simulations per parameter, for eight parameters). This analysis was performed for eight river basins (31 data points). The river basins used for the OAT analysis were the Pasig, Tullahan, Meycauayan, Jones Falls, Motagua, Ciliwung, Pesanggrahan, and Tiber (characteristics listed in table S6). We used correlation, regression analysis (coefficient of determination and Pearson correlation coefficient), and the ratio between the total simulated and observed plastic mass as evaluation metrics suggested by a study on sensitivity analysis workflows for environmental models (56). On the basis of the evaluation metrics, the sensitivity of the model was assessed, discussed, and lastly visualized in graphs and scatterplots.

Calculation process sensitivity analyses

For the AAT analysis, we performed 100 model runs for the three selected river basins, resulting in 100 output sets (plastic emission, ME) for three data points (rivers). This model output was compared with the model output of the best calibrated global parameter set for these three specific rivers. In the next step, to evaluate the performance and sensitivity of the complete parameter set, we calculated the coefficient of determination (r2) of the output of the 100 randomly generated parameter sets with the output of our best calibrated global parameter set (three points per set). Subsequently, we calculated the Pearson correlation coefficient between each parameter and the corresponding coefficient of determination (r2), and we calculated the Pearson correlation coefficient between each parameter and the corresponding plastic emission ME for each river basin. These two metrics provide insight into the impact of an individual parameter on the model output (plastic emission) and on the performance of the parameter variations to our best calibrated model on field observations, while other parameters are varied simultaneously. We divided the 100 parameter sets in three batches (R1 = 50 runs, R2 = 25 runs, and R3 = 25 runs) and assessed the impact of the variations of parameter values (input) on the model output for individual river basins using the method of visual inspections of the scatterplots (56).

For the OAT analysis, we performed 16 model runs for nine selected rivers with 31 data points (monthly and annual data, varying per river). Each parameter was increased or decreased by a factor of 10. The output (tons of plastic emission per river) was compared to both our best calibrated model and field observations. Subsequently, the coefficient of determination for the individual 16 scenarios and the best calibrated scenarios was calculated compared to the observations. Last, the ratio between the sum of modeled and observed plastic emissions was calculated for 16 scenarios and the best calibrated scenario.

Sensitivity analyses discussion

Figure S6 and table S17 show that the model is sensitive to all parameter values (adjusting parameter values leads to changes in model output) but not in an equal manner. Increasing probabilistic parameters results in increased model output; however, the r2 can decrease for both increasing and decreasing parameter values. Low correlation values (ranging between −0.10 and 0.10, for example) imply that variations (within bandwidths listed in table S15) have limited impact on the output of the model, while high correlation values (between −0.50 and −1.00 and between 0.50 and 1.00) imply that variations of parameter values have a greater impact on the model output and, therefore, on the coefficient of determination with the best calibrated scenario. After the first 50 runs (R1), we fixed the lower threshold stream order probability and lower threshold discharge probability from Eq. 5. These parameters show strong positive correlations with the r2, and hence, variations in these parameters have a large impact on the model output. This confirms the sensitivity of the model for river transport, P(O) in the framework, and supports assigning a realistic bandwidth to the parameters (as supported by expert elicitation). The next 25 runs (R2) reveal a high sensitivity for the stream order (and, to a lesser extent, discharge) coefficients, the coefficients that cause the threshold probability to increase for bigger rivers. Increasing the stream order coefficient results in a strong negative r2. Parameters are included in Eqs. 4 and 5, further described in table S9. On the basis of individual assessments of river basins, this parameter was fixed at a lower value, according to the best calibrated scenario, confirming both the sensitivity and the selected confined bandwidth. Last, R3 further confirms sensitivity to remaining parameters.

Responses for individual river basins can be different, for example, river basins with higher precipitation rates (Tullahan) have a higher correlation between the rainfall coefficient and model output. Scatterplots between the coefficient of determination of the model output and parameter variations for each parameter are presented in figs. S7 and S8. These scatterplots show that increasing parameter values results in higher model output for most parameters. This correlation depends on the parameter (sensitivity) and the river basin characteristics. The optimum value for a parameter can be assessed with the r2, showing that increasing parameter values can lead to a decreased performance. The longest river, the Motagua River, is very sensitive to the probability for plastic to be transported downstream.

The one-order-of-magnitude OAT sensitivity analysis shows that the model is sensitive to all parameter variations (fig. S9 and table S18); however, there are substantial differences between the r2 and ratio (modeled divided by observed total plastic emissions) per parameter. Increasing or decreasing the plastic mobilization probability P(M), by altering the rainfall and wind coefficient, reduces the r2 from 0.76 (best calibrated scenario) to between 0.28 and 0.66. The total plastic emission ratio varies between 0.23 and 2.40. These results suggest that increasing or decreasing wind and rainfall mobilization by a factor of 10 results in under- or overestimations of the model between a factor of 2.4 and 4. Variations of an order of magnitude in the probability for waste to reach a river P(R), by changing the slope and land use parameters, yield an r2 between 0.45 and 0.82 and a ratio of emission between 0.46 and 1.46. The model is most sensitive to changes in probability of waste in a river to reach the ocean, P(O). A decrease of an order of magnitude of the lower thresholds of stream order and discharge probability results in a negative r2 and large ratios between total observed and simulated plastic mass. The sensitivity of the model in correlation (r2) and the ratio between simulated and observed plastic emissions are presented per parameter and grouped for mobilization and for overland and riverine transport of plastic. Variations of parameters by a factor of 10 result in small ratios in output for mobilization and overland transport. For riverine transport, a decrease in the threshold for the probability of plastic to be transported 1 km downstream by a factor of 10 (from ~99 to 9.9%, for example) results in a factor of ~350 lower modeled transport. Calibration on field observations and parameter assessments from expert elicitation suggest river transport probabilities between 80 and 100%.

The sensitivity analyses show that the model output is sensitive to all parameters, both when they are varied simultaneously in interaction with other parameters (AAT) and individually (OAT). The AAT analysis indicates that the sensitivity also depends on river basin characteristics. The model performance decreases when parameter values move outside the bandwidth provided by expert elicitation and calibration, when compared to field observations. Both analyses suggest that the model is the most sensitive to the probability for plastic to be transported downstream a river.

Validation

We compared our best calibrated model simulation with independent field observations (tables S2 and S3). Data points were collected from published literature on (macro)plastic transport. These 84 data points originate from 51 rivers in six countries. We consider the Besos River (10 data points) an outlier because there are four weirs that may act as a sink for plastics (27) directly upstream of the observation point (28). The remaining 74 data points are roughly within one order of magnitude from modeled values, and we demonstrated a good order-of-magnitude correlation (coefficient of determination, r2 = 0.74, n = 74). The ratio between the sum of the modeled and observed data points is 1.73. Data points range from rivers with small emissions (~0.1 MT month–1) to large emissions (~10,000 MT month–1) (fig. S2). We do recognize that most of the data points used for validation are situated in Japan (39 of 51 rivers). This study reports combined microplastic and macroplastic emissions per prefecture with macroplastics representing, by far, most of the emitted mass. We conducted a separate assessment where we merged the 39 Japanese rivers into 1 data point, resulting in 36 data points from 12 rivers and 1 data point representing the Japanese prefectures (fig. S3). By doing so, the coefficient of determination became higher with r2 = 0.84, indicating that our validation results are not biased by the inclusion of many rivers from one single country.

Confidence intervals

We constructed confidence intervals for our model by considering the difference between 125 observed and simulated data points. We take the residuals of the logarithm of model simulations and field observations according toLog10(observation)Log10(model)=Log10(observation/model)(7)

We investigated whether the distribution of the ratio of the logarithm of 125 observations divided by model simulations can be assumed to be normally distributed. The Anderson (table S19) and Shapiro tests do not reject the normality assumption [Shapiro P(95) = 0.266 and Anderson-Darling test statistic = 0.259], and the Q-Q plot (fig. S10) shows a normal distribution. Therefore, the mean of the ratio of residuals and confidence intervals can be calculated with the SD for 68 and 95% (means ± 1 or 2 SDs, respectively). On the basis of 125 observations, the confidence interval for our model is a factor of 4 for the 68% confidence interval and 10 for the 95% confidence interval. It should be noted here that a negative ratio of −4 is equal to a ratio of 0.25 and implies that the model overestimates the plastic emission by a factor 4. Positive ratios mean that the model underestimates plastic emissions. Both a factor of −1 and 1 mean that observed and model-simulated plastic emissions are equal.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/7/18/eaaz5803/DC1

https://creativecommons.org/licenses/by-nc/4.0/

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank the anonymous peer reviewers for suggestions, which improved this manuscript. We thank the 60 participants of the expert elicitation. We thank I. Smith and R. Richardson for improving the English language of the manuscript. We are indebted to all people who collected field data on plastic in rivers and to the authors who published papers on plastic in rivers. We are grateful to nonauthoring contributors D. Higgins, R. Correia, K. van Oeveren, and M. Loozen for contribution to the sensitivity analyses and data collection and S. Pinson and R. de Vries for contribution to statistical analyses. We thank all the partners of The Ocean Cleanup involved in executing and enabling fieldwork and data collection. We thank B. Slat for fostering an environment that facilitated the freedom and resources to conduct the research required for this manuscript and for his contributions to the study itself. Funding: We would like to thank the supporters of The Ocean Cleanup who made the realization of this study possible. Author contributions: L.J.J.M., T.v.E., and R.v.d.E. designed the study. T.v.E. and L.J.J.M. conducted field expeditions to collect data. L.J.J.M. developed the model, and T.v.E., C.S., and L.L. reviewed the model. L.J.J.M., T.v.E., and L.L. wrote the manuscript. L.J.J.M. and L.L. prepared the figures. All authors reviewed the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: Data will be made available on Figshare and can be explored at www.theoceancleanup.com/1000rivers. All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article