Sampling

From July 27th to September 19th 2015, a total of 652 surface net tows were carried out between 25°N–41°N and 129–156°W by 18 participating vessels. In October 2016, we revisited our study area by conducting two flights with a Hercules C-130 aircraft that collected aerial imagery (n = 7,298 single-frame mosaics) to better quantify the larger and rarer >50 cm plastic objects (Fig. 1).

Figure 1 Field monitoring effort. Vessel (grey and dark blue lines) and aircraft (light blue lines) tracks and locations where data on buoyant ocean plastic concentrations were collected (circles). Grey circles (n = 350) represent areas sampled with a single Manta net tow by 17 participating vessels, between July and September 2015. Dark blue circles (n = 76) represent areas sampled with paired Manta and paired Mega net tows by RV Ocean Starr, between July and August 2015. Light blue circles (n = 31) show locations of RGB geo-referenced mosaics collected from a C-130 Hercules aircraft, in October 2016. This map was created using QGIS version 2.18.1 (www.qgis.org). Full size image

Vessels carried out net tows of 0.35–4 hours duration, while navigating at 0.7–6.8 knots. All trawls were designed to move away from the vessel to avoid wake effects on the capture efficiency of the devices. All vessel crews were trained with online material and one-to-one workshops that had been conducted prior to departure. While towing the trawl, the most experienced sailor aboard the vessel estimated the sea state (Beaufort scale) by measuring wind speeds and observing wave heights. This data was recorded in the standard datasheets provided, alongside the date, duration, as well as initial and final coordinates of each tow. The location and length of all net tows were confirmed during the post-processing phase by inspecting the position data from GPS trackers installed on all participating vessels. Most sampling stations encompassed a single net tow (n = 350 sampling stations) using a Manta trawl (0.5 mm square mesh, 90 cm × 15 cm mouth), which is one of the standard devices for quantifying plastic pollution levels. With the largest participating vessel (RV Ocean Starr), we simultaneously towed two Manta trawls, alongside two large Neuston trawls (1.5 cm square mesh, 6 m × 1.5 m mouth, of which 0.5 m above the water line; thereafter called ‘Mega trawls’) at every sampling location (n = 76 stations). After each Manta net tow, the net was rinsed from the outside with seawater, and its single-use cod-end removed, closed with staples and placed in an individual zip-lock bag. After each Mega trawl tow, the net was also rinsed from the outside with seawater and its large cod-end opened in a box filled with seawater. All buoyant plastics were then removed, wrapped in aluminium and placed in labelled plastic bags. The whole content captured by the Manta trawls was stored, while the organisms captured by the Mega trawls (mostly alive) were released back into the ocean. All samples were stored in a fridge or freezer while at-sea, and in a FedEx cool box (2–8 °C) or reefer (−2 °C) while being shipped to the laboratory. Even though we were careful when handling samples, some debris items may have been broken during transportation, leading to some bias in our debris size distribution. Detailed information related to these net tows (i.e. coordinates, metocean conditions, sampling times and durations) is provided in Figshare33.

The aerial surveys sampled a far greater area (311.0 km2) than the trawl surveys described above (3.9 km2 and 13.6 km2, for Manta and Mega net tows, respectively), thus yielding a more reliable quantification of debris larger than 50 cm, which are relatively rare. Both flights started and ended at Moffett Airfield near Mountain View, California. The first aerial survey was conducted on October 2nd 2016 sampling from 18:56 to 21:14 UTC time, at a constant latitude of 33.5°N, and longitudes varying from 141.4°W to 134.9°W. The second survey started on October 6th 2016 sampling from 22:14 to 0:37 UTC, from 30.1°N, 143.7°W to 32.9°N, 138.1°W. While in survey mode, the aircraft flew at an altitude of approximately 400 m and at a ground speed of 140 knots. Sampling transects targeted areas where the sea state conditions were the lowest, based on the weather forecast, including sea surface atmospheric pressure, cloud cover, wind speed at 10 m above sea level and boundary surface layer height provided by NOAA’s Global Forecasting System, as well as significant wave height and peak period data distributed by NOAA’s WaveWatch3 model outputs. Even though we surveyed floating debris using trained observers and three types of sensor (Lidar, SWIR imager, and RGB camera), here we only analyse information coming from the geo-referenced mosaics produced by a RGB camera (CS-4800i) that generally took photographs every second during surveying time (frame size = ~360 m across track, ~240 m along track, ~0.1 m resolution).

Trawl samples processing

Trawl samples were separately washed into a sieve tower (five Glenammer Engineering Ltd sieves, with 0.05 cm, 0.15 cm, 0.5 cm, 1.5 cm, and 5 cm square apertures) that divided the material into the following size classes: 0.05–0.15 cm, 0.15–0.5 cm, 0.5–1.5 cm, 1.5–5 cm, and >5 cm. Debris items >5 cm were then manually sorted into 5–10 cm, 10–50 cm, and >50 cm classes by measuring the object lengths (widest dimension of the object) with a ruler. Buoyant debris was separated from biomass by placing the material within each sieve in filtered saltwater (salinity 3.5%, temperature 19–23 °C). Lab personnel stirred the material many times to insure floating particles were detached from the biomass material. Floating objects identified as buoyant debris were manually extracted from the water surface using forceps, separated into types, and counted. Buoyant debris was classified into material type (plastic, glass, paraffin, tar, rubber, wood, pumice, seed or unknown), with plastics being further divided into the following categories: (1) ‘H’ type – fragments and objects made of hard plastic, plastic sheet or film; (2) ‘N’ type – plastic lines, ropes, and fishing nets; (3) ‘P’ type – pre-production plastic pellets in the shape of a cylinder, disk or sphere; and (4) ‘F’ type – fragments or objects made of foamed material (e.g. expanded polystyrene). Once counted and categorized, the pieces were washed with distilled water, transferred to aluminium dishes, dried overnight at 60 °C, and weighed using an OHAUS Explorer EX324M (0.0001 g readability) for objects <5 cm, and a OHAUS Explorer EX12001M (0.1 g readability) for objects >5 cm.

To best characterize the ocean plastic accumulating within the GPGP, we performed additional analyses with the material collected. Firstly, 10 pieces within each plastic size/type category (n = 220 pieces) were selected for polymer composition analysis by Fourier-transform infrared spectroscopy (FT-IR). The readings were done using a Perkin Elmer Spectrum 100 FT-IR equipped with a universal ATR accessory (range = 600–4000 cm−1). The respective polymer type was determined by comparing sample FT-IR spectra against known spectra from a database (Perkin-Elmer ATR of Polymers Library). Secondly, we screened all plastic debris collected for production dates, as well as any writings statements giving information on its origin (i.e. language and ‘made in’ statements). Lastly, we classified plastic items from ‘H’ and ‘L’ types collected at 30 RV Ocean Starr stations into object types (e.g. bottle lids, bags, bottles, etc). As ‘H’ objects larger than 50 cm were relatively rare, we analysed 10 extra RV Ocean Starr stations for this type/size category. If the object type of a fragment could not be determined, we classified the piece as either hard plastic fragment or film fragment depending on its wall thickness and flexibility34. We used Manta trawl samples to characterise objects within size classes 0.15–0.5 cm, 0.5–1.5 cm, and 1.5–5 cm, and Mega trawl samples to characterise objects within size classes 5–10 cm, 10–50 cm, and >50 cm. Plastics within our smallest size class (0.05–0.15 cm) were not considered in this ‘Object Type’ analysis due to the difficulty of handling and identifying small fragments.

The numerical/mass concentrations of buoyant plastic items (count/kg of plastic per km2 of sea surface) measured by each net tow were calculated for all plastic size/type categories separately. To do so, we divided the count and weight of plastic objects within each category by the towed area of the sample. We calculated the towed area by multiplying net mouth width (90 cm for Manta trawl, 6 m for Mega trawl) by tow length (determined from GPS position data). The average area covered by Manta net tows was 0.008 km2 (SD = 0.004, min–max: 0.001–0.018 km2), while the average area covered by Mega net tows was 0.090 km2, (SD = 0.013, min–max: 0.046–0.125 km2). As buoyant plastics can be missed by surface trawls due to wind-driven mixing, we then estimated the ‘depth-integrated’ mass and numerical plastic concentrations (Ci) for all type/size categories at each of the trawl sampling locations using the equations described in ref.35. Supplementary Methods 1 provides details on how Ci was calculated as a function of ocean plastic terminal rising velocity (Wb), depth sampled by the trawl, and sea state. It also describes how we measured Wb for each of the type/size categories of this study. After comparing plastic concentration results obtained by paired Manta and Mega net tows (n = 76 locations), we decided to use Manta and Mega trawl samples to quantify debris 0.05–5 cm and 5–50 cm in size, respectively. The comparison results and reasoning behind such decision in provided in Supplementary Methods 2.

Aerial imagery processing

All RGB images taken during our survey flights (n = 7,298) were georeferenced using accurate aircraft position and altitude data collected during the surveys. They were then inspected by two trained observers and a detection algorithm. Observers inspected all images at full-screen on a Samsung HD monitor (LU28E590DS/XY) and those single-frame mosaics containing debris were uploaded into QGIS software (Version 2.18.3–Las Palmas) to record their position and characteristics. We trust we had a very small number of false positives and a high number of false negatives. This is because the observers took a conservative approach: they only logged features as debris when they were very confident with its identification. As such, many features that could be debris, but resembled other natural features, such as sun glint and breaking wave, were not logged into our ocean plastic dataset. Once this work was finalised, we ran an experimental algorithm capable of detecting potential debris in all our RGB mosaics as a quality control step. To avoid any false positives, all features detected by the algorithm were also visually inspected by an observer and only those visually identified as debris were logged in our QGIS database. For every sighting, we recorded position (latitude, longitude), length (widest dimension of the object), width, and object type: (1) ‘bundled net’ – a group of fishing nets bundled tightly together; they are commonly colourful and of a rounded shape; (2) ‘loose net’ – a single fishing net; they were generally quite translucent and rectangular in shape; (3) ‘container’ – rectangular and bright objects, such as fishing crates and drums; (4) ‘rope’ – long cylindrical objects around 15 cm thick; (5) ‘buoy/lid’, rounded bright objects that could be either a large lid or a buoy; (6) ‘unknown’ – objects that are clearly debris but whose object type was not identified, they were mostly irregular-shaped items resembling plastic fragments; and (7) Other – only one object was successfully identified but did not belong to any category above: a life ring. We recorded 1,595 debris items (403 and 1,192 in flights 1 and 2 respectively); 626 were 10–50 cm and 969 were >50 cm in length. Most of them were classified as ‘unknown’ (78% for 10–50 cm, 32% for >50 cm), followed by ‘buoy or lid’ (20%) and ‘bundled net’ (1%) for 10–50 cm debris, and by ‘bundled net’ (29%), ‘container’ (18%), ‘buoy or lid’ (9%), ‘rope’ (6%), and ‘lose net’ (4%) for >50 cm debris. To calculate ocean plastic concentrations, we grouped the geo-referenced images into 31 ~10 km2 mosaics. For numerical concentrations, we simply divided the number of debris pieces 10–50 cm and >50 cm within each mosaic by the area covered. To estimate mass concentrations, we had to first estimate the mass of each object spotted, then we separately summed the mass of 10–50 cm and >50 cm debris within each mosaic by the area covered. More information on how we estimated the mass of each spotted objects is provided in Supplementary Methods 3.

Numerical model formulation

Ocean plastic pathways can be represented by Lagrangian particle trajectories31. In our framework, particles were advected by the following environmental drivers: sea surface currents, wave induced Stokes drift and winds. Starting from identical particle releases, we produced a series of forcing scenarios to represent the diversity in shape and composition of ocean plastics. Starting from using sea surface current only, we gradually added forcing terms representing the actions of atmospheric drag and wind waves on buoyant debris. The action of wind was simulated by considering the displacement of particles as a fraction of wind speed at 10 m above sea level. This is referred as the ‘windage coefficient’. We assessed different windage coefficient scenarios including 0%, 0.1%, 0.5%, 1%, 2% and 3%. We sourced global sea surface currents (1993 to 2012) from the HYCOM + NCODA global 1/12° reanalysis (experiment 19.0 and 19.136,37,38), and wind (10 m above sea level) speed and direction data (1993 to 2012) from NCEP/NCAR global reanalysis39. Wave induced Stokes drift amplitude was calculated using wave spectrum bulk coefficients (significant wave height, peak wave period and direction) from Wavewatch3 model outputs40.

For every forcing scenario, particles were identically and continuously released in time from 1993 to 2012 following spatial distributions and amplitudes of significant ocean plastic sources on land (coastal population hotspots23 and major rivers24) as well as at sea (fishing26,41, aquaculture42 and shipping industries43). Source scenarios were combined using relative source contribution as well as geographical distribution presented in Supplementary Methods 4. We advected global particles in time using the forcing scenarios described above and successfully reproduced the formation of oceanic garbage patches, with the shape and gradient of particle concentrations in these areas differing amongst forcing scenarios. We computed daily particle visits over 0.2° resolution grids corresponding to our observation domain and extending from 160°W to 120°W in longitude and 20°N to 45°N in latitude. The number of daily particle visits was uniformized over the total number of particles present in the global model at a given time. The model-predicted non-dimensional concentration δ i of cell i, was calculated as follows:

$${\delta }_{i}=\sum _{s}{\alpha }_{s}{\delta }_{i,s}$$ (1)

where α s is the non-dimensional weight relative to the contribution of source s and δ i,s is the percentage of global particles from source s in cell i. δ i,s is calculated with the number of particles n i,s from source s in cell i over the total number of global particles Σ i n s from source s:

$${\delta }_{i,s}=\frac{{n}_{i,s}}{{\sum }_{i}{n}_{s}}$$ (2)

Numerical model calibration

We collected measurements at sea in 2015 and 2016, but our numerical model uses ocean circulation reanalysis covering the period from 1993 to 2012. Modelled ocean circulation data post-2012 is available from HYCOM however not as a reanalysis product. As such, we decided not to use it in this study. As initial model particles released in 1993 significantly start to accumulate in the area after about 7 years, we averaged uniformized daily particle visits over 12 years, from 2000 to 2012. We grouped observed debris size classes in four categories: microplastics (0.05–0.5 cm), mesoplastics (0.5–5 cm), macroplastics (5–50 cm) and megaplastics (>50 cm). We compared model predictions against depth-integrated microplastic concentrations as this dataset collected by Manta trawls had the largest spatial coverage. Mass concentrations derived from trawl measurements were grouped in 0.2-degree resolution cells and compared against model-predicted non-dimensional concentration δ for the five different forcing scenarios. The best model fit was found for the forcing scenario with sea-surface current only (R2 = 0.52, n = 277 cells). The regression coefficient declined as we increased the atmospheric drag term (R2 = 0.39 to 0.21 depending on windage coefficient).

As we analysed the accumulation of model particles in the GPGP region, we noticed significant seasonal and inter-annual variations of the GPGP position. The modelled GPGP dimensions was relatively consistent throughout our 12 years of analysis, but the relative position of this accumulation zone varied with years and seasons. We first decided to test our model for seasonal variation by comparing our microplastic concentrations (measured in July – September 2015) against modelled concentrations averages for the July–September periods of 2000 to 2012. This comparison yielded poorer results (R2 = 0.46 to 0.21, depending on forcing scenario) than with the 12 years average solution (R2 = 0.52) as the July–September GPGP position varied substantially among years.

The relation between the accumulation of marine debris in the North Pacific and climate events such as El Niño Southern Oscillation (ENSO) and the Pacific Decadal Oscillation (PDO) has previously been discussed18. As such, to account for inter-annual variation, we compared latitudinal and longitudinal position of the GPGP against these two climate indexes: ENSO and PDO. We found that 2002 and 2004 were similar to the conditions experienced during our multi-vessel expedition. Thus, we compared our measurements against particle visit averages for July–September of 2002 and 2004 combined. This second attempt exhibited better results (R2 = 0.58 to 0.41, depending on forcing scenario), suggesting that climatic events such as ENSO or PDO influence the average position of the GPGP. Therefore, we decided to use the July–September average for 2002 and 2004, which better accounts for inter-annual variations in the GPGP position. More information on the selection of years for calibrating the model against trawl and aerial surveys data is provided in Supplementary Methods 5. The best fit between model predictions and microplastic observations was found once again for the forcing scenario with sea surface current only (R2 = 0.58, n = 277). The best regression fit between measured and modelled microplastic concentrations had a = −8.3068 and b = 0.6770 in the parametric formulation:

$$\,{c}_{mod}=\,{10}^{\frac{{\mathrm{log}}_{10}\delta -a}{b}}$$ (3)

From this formulation, we computed the modelled microplastics mass concentration in our domain area and extracted contour levels by order of magnitude, from 0.01 g km−2 to 10 kg km−2. The GPGP as defined in this study corresponds to the 1 kg km−2 microplastic mass concentration level covering an area of 1.6 million km2 and depicted as a bold line in Fig. 2a. As a validation, we categorized microplastics measurements inside and outside the 1 kg km−2 contour line (Fig. 2b). For stations inside the model-predicted GPGP, the median measured microplastic concentration was 1.8 kg km−2 (25th–75th percentiles = 3.5–0.9 kg km−2) while for stations outside, the median was 0.3 kg km−2 (25th–75th percentiles = 0.2–0.7 kg km−2). Using our calibrated microplastics distribution, we computed the mass and numerical concentration for individual size classes from scaling modelled concentrations by the ratio between average modelled microplastics distribution inside the GPGP and averaged measured concentrations per size class of stations inside the patch. A comparison between measured and modelled mass/numerical concentrations for all ocean plastic size classes is given in Fig. 2c and d.

Figure 2 Numerical model calibration. (a) The GPGP boundary (blue line) is estimated by comparing microplastic concentration measurements (circles) to model particle visit averages that accounted for seasonal and inter-annual variations. This map was created using QGIS version 2.18.1 (www.qgis.org). (b) Model validation showing median measured mass concentration for microplastics of stations outside and inside our predicted 1 kg km−2 GPGP boundary. Bars extend from 25th to 75th percentile while whiskers extend to minimum and maximum non-outlier. Outliers are represented as crosses. (c) Measured mass concentrations versus modelled mass concentrations for microplastics, mesoplastics, macroplastics and megaplastics. (d) Same as (c) but with numerical concentrations. Full size image

Our confidence intervals were formulated to account for uncertainties in both sampling and modelling. For the trawl collection (i.e. micro-, meso- and macroplastics), we considered uncertainties related to the vertical-mixing corrections applied to surface concentrations using reported sea state and plastic’s rising velocities (see Supplementary Methods 1). For the aerial mosaics, we accounted for uncertainties related to estimating the mass of sighted objects based on correlations between top-view area and dry weight of objects collected in the trawls (see Supplementary Methods 3). Finally, to account for modelling uncertainties, we added (respectively subtracted) the standard error of measured concentration to (resp. from) the mean upper (resp. lower) mass concentration when scaling the microplastics distribution to individual size classes.

Characterisation by types, sources and forcing scenarios

The total estimated mass load of ocean plastic in the GPGP by size classes were further divided by types. We calculated the average fraction in mass of individual ocean plastic types per sampling event for stations inside the patch (Supplementary Table 1) and derived the contribution of types ‘H’, ‘N’, ‘F’ and ‘P’. Further, as we predominantly observed debris originated from marine sources, we investigated the source contribution predicted by our calibrated modelled distribution. For individual model cells, we calculated the percentage of Lagrangian particle visits from individual sources. As initial particles were weighted in accordance to estimated global inputs, model particles from marine sources originally represented 28.1% of the total amount of material with fishing (17.9%), aquaculture (1.3%) and shipping (8.9%). We calculated the difference from this initial percentage value for each model cell and reported it to the predicted total mass concentration. In doing so, we defined ‘anomalies’ in marine source contribution in the North Pacific and expressed these in unit of mass per surface area. Finally, even though our calibrated model considered sea surface current only, we compared the predominance of forcing scenarios by evaluating the respective number of particle visits for each model cells. We computed contours around the GPGP for individual forcing scenarios in a way that the material contained inside each contour is equal to our initial forcing scenario (i.e. sea surface current only).

The dependence of the particle trajectory on the windage coefficient predicted by our model is in good agreement with sightings and modelling of debris originated from the 2011 Tohoku tsunami in Japan44,45. The first identified Japanese debris items that arrived after 10 to 12 months on the North American shores were objects with high windage such as buoys, boats and floating docks. Debris also arrived on Hawaiian Islands 18 months after the incident. The time of arrival was closely related to object types, starting in the first year with large oyster farm buoys and other floats, containers, and canisters. In the second year more buoys, tipped boats, fridges, and pallets arrived, followed later by timber beams and wooden debris. Our model predicted that only objects with a windage coefficient above 3% could arrive on Hawaii in the second year after the 2011 tsunami. Objects with windage coefficient ranging from 1 to 2% would reach Hawaii during the third year, while objects with no windage would mostly accumulate in the GPGP, north east of the archipelago.

Long term analysis

The definition of a dynamic GPGP boundary that accounts for seasonal and inter-annual variabilities allowed us to estimate which sea surface trawl data points from the literature are inside or outside the GPGP region. Therefore, we used our calibrated model to assess the decadal evolution of microplastic mass concentrations (kg km−2) within and around the GPGP. Concentration data from the literature (Supplementary Table 2) was obtained from published datasets or digitized from figures when not available digitally17,46,47. When data was reported in unit of mass per volume of water48, we used the net tow depth to calculate the concentration per surface area unit. When only numerical concentration was reported22,48, we estimated mass concentration by using the average ocean plastic mass from net tows where both mass and numerical concentrations were reported (m = 3.53 mg, SE: 0.10 mg, n = 872).

We compared the model-predicted GPGP boundary with the locations of samples collected between 1999 and 201221,22,48,49. Samples collected before 199917,46,47,48 were compared against the GPGP position estimated for the sampled months and years in the 1999–2012 period that had similar ENSO and PDO values (See Supplementary Methods 6). Using our dynamic GPGP model boundary as reference, we classified each net tow into 3 categories: (1) sampled within the GPGP boundary, (2) sampled outside the GPGP boundary, but above 20°N and below 45°N and (3) sampled in the rest of the North Pacific. We only used net tows from the first two categories above so that concentration statistics for outside the patch were not biased by measurements taken in equatorial and polar waters, where concentrations were very low. We then grouped these microplastic concentration observations from plankton net trawls by decades, taking data recorded between 1965–1974 (n = 20 inside and n = 58 outside17,48), 1975–1984 (n = 0 inside and n = 19 outside46), 1985–1994 (n = 4 inside and n = 2 outside47), 1995–2004 (n = 2 inside and n = 252 outside22,49), 2005–2014 (n = 195 inside and n = 861 outside21,22,48) and finally 2015 (n = 288 inside and n = 213 outside; this study). We calculated the mean (± standard error) of measured microplastic mass concentration per decades for within and around the GPGP boundary. Finally, we extracted decadal trends by fitting an exponential function (R2 = 0.94) assuming null concentrations at the beginning of the 20th century. The exponential fit exhibited better results than linear, quadratic or cubic functions (R2 = 0.71, R2 = 0.86 and R2 = 0.91, respectively).