Modelling floating riverine litter in the south-eastern Bay of Biscay: a regional distribution from a seasonal perspective

Although rivers contribute to the flux of litter to the marine environment, estimates of riverine litter amounts and detailed studies on floating riverine litter behaviour once it has reached the sea are still scarce. This paper provides an analysis of the seasonal behaviour of floating marine litter released by rivers within the south-eastern Bay of Biscay based on riverine litter characterizations, drifters, and high-frequency radar observations and Lagrangian simulations. Virtual particles were released in the coastal area as a proxy of the floating fraction of riverine litter entering from rivers and reaching the open waters. Particles were parameterized with a wind drag coefficient (C_d) to represent their trajectories and fate according to the buoyancy of the litter items. They were forced with numerical winds and measured currents provided by high-frequency radars covering selected seasonal week-long periods between 2009 and 2021. To gain a better insight into the type and buoyancy of the items, samples collected from a barrier placed at the Deba River (Spain) were characterized at the laboratory. Items were grouped into two categories: low-buoyancy items (objects not exposed to wind forcing, e.g. plastic bags) and highly buoyant items (objects highly exposed to wind forcing, e.g. bottles). Overall, low-buoyancy items encompassed almost 90 % by number and 68 % by weight. Weakly buoyant items were parameterized with C_d = 0 % and highly buoyant items with C_d = 4 %; this latter value is the result of the joint analysis of modelled and observed trajectories of four satellite drifting buoys released at the Adour (France), Deba (Spain), and Oria (Spain) river mouths. Particles parameterized with C_d = 4 % drifted faster towards the coast through the wind, notably during the first 24 h. In summer, over 97 % of particles beached after 1 week of simulation. In autumn this value fell to 54 %. In contrast, low-buoyancy items took longer to arrive at the shoreline, particularly during spring with fewer than 25 % of particles beached by the end of the simulations. The highest concentrations (>200 particles km⁻¹) were recorded during summer for C_d = 4 % in the French region of Pyrénées-Atlantiques. Results showed that the regions in the study area were highly affected by rivers within or nearby the region itself. These results couple observations and a river-by-river modelling approach and can assist decision-makers on setting emergency responses to high fluxes of floating riverine litter and on defining future monitoring strategies for heavily polluted regions within the south-eastern Bay of Biscay.

1 Introduction

Rivers act as key vectors bringing improperly disposed and mismanaged litter from land into the marine environment. Riverine litter poses a large threat to freshwater systems by degrading aquatic life, impacting freshwater quality, and increasing economic losses linked with human activities (van Emmerik and Schwarz, 2020; Al-Zawaidah et al., 2021). However, most of the litter research conducted to date has focused on marine environments (87 %) when compared to freshwaters systems (13 %) (Blettler et al., 2018). Indeed, riverine litter contributions to oceans are still uncertain, and results vary depending on the input data and the model applied (Lebreton et al., 2017; Schmidt et al., 2017; Mai et al., 2020). Recent findings derived from extensive modelling efforts suggest that about 1600 rivers worldwide account for 80 % of plastic inputs to the ocean, with small urban rivers among the most polluting (Meijer et al., 2021). Models require comprehensive field data and consistent and harmonized protocols to validate the amounts, type, and size of riverine inputs (González-Fernández and Hanke, 2017; Wendt-Potthoff et al., 2020; Margenat et al., 2021). Such comprehensive data were obtained in Europe thanks to the RIMMEL project (González-Fernández and Hanke, 2017). This research concluded that between 307 and 925 million floating riverine litter items are annually transferred into the ocean mainly through small rivers, streams, and coastal run-off (González-Fernández et al., 2021).

Once it has reached the sea, floating riverine litter can accumulate close to the shoreline or it can be transported to open waters, reaching even remote areas far from the coast. Indeed, the distribution and fate of floating litter in the marine environment is affected by the metocean conditions (currents, turbulence, wind) but also by the buoyancy of the objects (Ryan, 2015; Lebreton et al., 2019; Maclean et al., 2021). Objects with low buoyancy are mainly driven by currents, by contrast with highly buoyant items, which are driven along the water surface partially by winds. This wind effect (“windage”) on floating marine litter behaviour has been further investigated by Lagrangian modelling studies in the open ocean (Allshouse et al., 2017; Maximenko et al., 2018; Lebreton et al., 2019; Abascal et al., 2009) when compared to the coastal area (Critchell and Lambrechts, 2016; van Utenhove, 2019; Tong et al., 2021). The lack of observational data is one of the key limitations for parameterizing the windage effect and accurately predict floating marine litter behaviour. However, observations derived from drifting buoys, such as those provided for decades by the Global Drifter programme, have been used to fill this gap. They have allowed simulating more realistic floating marine litter pathways from origin to fate by integrating experimental windage parameterizations and the corresponding comparison between observed and modelled trajectories (Duhec et al., 2015; Pereiro et al., 2018; Rizal et al., 2021). Nowadays, more affordable and environmentally friendly solutions are gaining force among researchers, as drifters are built using biopolymers (Novelli et al., 2017; D'Asaro et al., 2020) or have compact and lightweight designs with a GPS-tracking component for easy deployment (Meyerjürgens et al., 2019; van Sebille et al., 2021).

At the coastal scale, windage parameterization combined with realistic knowledge on coastal circulation becomes crucial to reduce the uncertainties of modelled trajectories (Van Sebille et al., 2020). Land-based high-frequency radar systems (hereafter HF radars, Rubio et al., 2017) offer the opportunity to monitor surface currents in coastal areas, where the transport processes are significantly more complex than in open-ocean waters due to the effect of the coast, the bathymetry, and other local forcings (e.g. river discharges or coastal upwellings). In the south-eastern Bay of Biscay (hereafter SE Bay of Biscay), as part of the operational oceanography system EuskOOS, an HF radar provides near-real-time surface currents fields. The system has already been used to study surface coastal transport processes in the area in combination with multisource data (Manso-Narvarte et al., 2018, 2021; Rubio et al., 2011, 2013, 2018, 2020; Solabarrieta et al., 2014, 2015, 2016). The HF radar is also a good example of effective monitoring of surface currents with strong potential for floating marine litter management. Research conducted by Declerck et al. (2019) in the SE Bay of Biscay provided the first assessment of floating marine litter transport and distribution in the region, coupling surface current observations from the EuskOOS system, Lagrangian modelling, and riverine inputs. Nowadays, these observations are used by local authorities both in real time and in hindcast in the framework of the operational service FML-TRACK to collect floating marine litter in the area. However, the accurate modelling of the transport and fate of floating marine litter needs to consider the variety of floating objects and sources and additional physical parameterization, such as windage.

This paper aims at estimating the seasonal behaviour of the floating marine litter fraction released by rivers within the SE Bay of Biscay reaching open waters. To do so, a Lagrangian model was forced by real observations from the EuskOOS HF radar and particles were parameterized to represent floating marine litter trajectories of two groups of items according to their buoyancy. Riverine litter collected from a local barrier was characterized at the laboratory to explore the fraction of high- and low-buoyancy items. Since most of the items were low-buoyancy, simulations considering only surface currents were performed as the reference. Complementary Lagrangian simulations for highly buoyant items (and less abundant in the area) were also performed. In this case, four low-cost buoys with a similar buoyancy to certain highly buoyant items were built and released in three different rivers. Drifter data were used to parameterize the wind effect on this type of item and consequently achieve more accurate results.

2 Study area

The study was conducted in the SE Bay of Biscay, between north-eastern Spain (Basque Country) and south-western France (Landes). The study area extends from 43.27 to 44.58^∘ N and from 3.18 to 1.27^∘ W, falling within the coverage area of the HF radar station of the operational oceanography system EuskOOS (Fig. 1). The study area comprises two Basque regions (Bizkaia (Spain) and Gipuzkoa (Spain)), two French departments (Pyrénées-Atlantiques and Landes), and eight rivers (Deba (Spain), Urola (Spain), Oria (Spain), Urumea (Spain), Oiartzun (Spain), Bidasoa (Spain), Nivelle (France), and Adour (France)). The mean annual river discharge varies widely between rivers from 3.71 m³ s⁻¹ (Oiartzun) to 350 m³ s⁻¹ (Adour) (Sheppard, 2018), and the population density differs between the Spanish and French border: 44.8 inhabitants km⁻² (Landes) to 303.7 inhabitants km⁻² (Basque Country) (Eurostat, 2019). The bathymetry in the SE Bay of Biscay is characterized by the presence of a narrow continental shelf ranging 7 and 24 km wide in the Basque area, gradually increasing along the French coast up to about 70 km (Bourillet et al., 2006; Rodríguez et al., 2021). The continental shelf in the SE Bay of Biscay comprises two main areas: the Aquitaine shelf with a N–S orientation and the Cantabrian shelf with an E–W orientation. The continental slope is very pronounced, with a slope of the order of 10 %–12 % (Sheppard, 2018). Over the continental shelf, the ocean circulation is marked by seasonal variability. At shorter temporal scales, circulation in the study area is mostly modulated by the bathymetry and the coastal orientation, the density-driven currents, and winds (Le Boyer et al., 2013; Solabarrieta et al., 2014). Tidal currents are quite weakly constrained by the topography and the width of the continental shelf (Lavin et al., 2006; González et al., 2007; Karagiorgos et al., 2020). Along-shelf currents are more intense and persistent during winter and autumn (about 10–15 cm s⁻¹), contrary to the other seasons, especially in summer (about 2.5 cm s⁻¹) (Charria et al., 2013). In winter, the prevailing SW winds causes an E to N flow from the Spanish coasts towards the French coasts. The moderate to strong NW winds occurring in spring and summer induce a S and SW surface current circulation accompanied by a greater variability (Solabarrieta et al., 2015). In winter, westerly winds in the Basque coast reinforce the slope current (named the Iberian Poleward Current (IPC)), a warm and saline intrusion trapped within 50 km of the shelf edge, reaching its greatest velocities (up to 70 cm s⁻¹) during this season. The IPC favours the along-slope transport of water masses (Solabarrieta et al., 2014; Porter et al., 2016). The exchange between shelf and deep-sea waters in winter is associated with the generation of eddies, from the interaction of currents with the topography (Lavin et al., 2006; Rubio et al., 2018; Teles-Machado et al., 2016). Maximum run-offs combined with SW winds also allow river plumes to spread northwards and along the French shore during winter. However, this path changes in spring, when river discharges are reduced and winds blow from the north-west (Lavin et al., 2006; Puillat et al., 2006).

First global modelling studies coupling ocean circulation and Lagrangian particle-tracking models reported that the SE Bay of Biscay is a hotspot for floating marine litter (Lebreton et al., 2012; van Sebille et al., 2012). Recent Lagrangian modelling studies combining measured and predicted surface currents by the HF radar and the Iberian Biscay Irish System (IBI) Copernicus model revealed that floating marine litter circulation in the SE Bay of Biscay is marked by a high seasonal variability. Results showed a higher retention during spring and summer and a northward dispersion along the French coast during autumn and winter (Declerck et al., 2019; Rubio et al., 2020). Surface currents derived from the Regional Ocean Modelling System (ROMS) and a particle-tracking model were combined by Pereiro et al. (2019) to track the numerical drifters representing floating marine litter in the Bay of Biscay. In this study, longer residence times and higher concentrations were observed in the SE Bay of Biscay when compared to north-western Iberian coastal waters, particularly in winter. From numerical simulations run using the HYCOM model, Rodríguez-Díaz et al. (2020) showed that floating marine litter items with high windage (C_d = 3 %–5 %) tend to accumulate in nearshore areas of the Bay of Biscay or end up beached. This trend is consistent with recent numerical simulations combining surface currents from the operational IBI and the numerical model TESEO that also revealed that the highly buoyant items (C_d = 4 %) rapidly beach in the SE Bay of Biscay, mainly in spring and summer (Ruiz et al., 2022a). Since June 2020, innovative detection and tracking solutions combining ocean modelling and remote observation systems have been operating in the SE Bay of Biscay to support floating marine litter reduction strategies both downstream (interception at the sea with collection vessels and on beaches with cleaning facilities) and upstream (source identification and reduction) (Delpey et al., 2021). However, research on floating marine litter behaviour in the SE Bay of Biscay is still in its early stage. Further experiments are needed to fully understand the role of windage, waves, and tides in the complex 3D circulation patterns governing coastal accumulation.

Figure 1 Study area with the release locations of the satellite drifting buoys and the riverine barrier. Dots in light yellow represent the nodes of the HF radar grid. Dots in orange represent the trajectories of the buoys. Numbers with stars in pink correspond to the particle releasing location for floating marine litter simulations: (1) Deba, (2) Urola, (3) Oria, (4) Urumea, (5) Oiartzun, (6) Bidasoa, (7) Nivelle, and (8) Adour rivers.

3 Methods and data

3.1 Riverine litter sampling

In spring 2018, a riverine barrier was placed in the Deba River (Gipuzkoa) to retain and collect floating riverine litter during low to moderate flows. This barrier enabled a passive sampling to characterize litter items in the lab. The barrier, which consisted of a nylon artisanal net supported by hard floats (buoys), was 40 m long and 0.6 m high with a 60 mm mesh size (see photos in Appendix A). The sampling was conducted weekly from April to June 2018. In total eight riverine litter samples were collected. Litter items were quantified, weighed, and categorized in the lab according to the master list included in the Guidance on Monitoring of Marine Litter in European Seas (Joint Research Centre, 2014). Items were grouped into seven types of material (artificial polymer materials, rubber, cloth/textile, processed/worked wood, paper/cardboard, metal, and glass/ceramics) and further classified into 44 categories (see the classification in Appendix B). Riverine litter items were also categorized into two groups (low- and high-buoyancy items) considering their exposure to wind based on Ruiz et al. (2022a).

3.2 Drifter observations

Four satellite drifting buoys (herein after “low-cost buoys”) were built by the authors and deployed one by one in the river mouths of the Deba (Buoy A), Oria (Buoy B), and Adour (buoys C and D) between April 2018 and November 2018 (Fig. 1, Table 1). The low-cost buoys provided positioning every 5 min using satellite technology. Low-cost buoys were 9 cm in height and 9.5 cm in float diameter and weighed approximately 200 g (Fig. 2). A GPS (SPOT Trace device) powered by four AAA cells was placed in the bottom of a high-density polyethylene (HDPE) plastic container sealed to guarantee water tightness. They were chosen because of their capability to ensure a reasonable balance between an accurate signal emission and their purchase and communication fees. SPOT Trace devices have been used over the past few years in coastal and open-ocean applications in a wide range of studies. Studies range from calibrating HF radars (Martínez Fernández et al., 2021) and tracking drifting objects such as icebergs (Carlson et al., 2020), pelagic Sargassum (Putman et al., 2020; van Sebille et al., 2021), or fishing vessels (Widyatmoko et al., 2021; Hoenner et al., 2022) to search and rescue training (Russell, 2017) and oil spill and litter monitoring (Novelli et al., 2018; Meyerjürgens et al., 2019). Almost two-thirds of the buoy floated above the water surface, thus preventing any satellite signal losses. Buoys A and D and transmitted their positions on an ongoing basis until their landing. Buoys B and C stopped emitting while they were drifting. In all cases, battery lifetime was enough for an adequate performance of the buoys. Once on land, citizens collected the buoys and reported their corresponding location.

Figure 2 Main components of the low-cost buoy. The structure: (a) HDPE container and SPOT Trace device powered by four AAA cells. Assembly process: (b) final appearance once the buoy is sealed; the buoy is labelled with contact information both within and outside; (c) the SPOT Trace was fixed at the base of the container with adhesive tape to avoid twists and turns of the buoy.

Table 1 Locations, periods, and distances covered by the drifting buoys.

3.3 HF radar current observations and wind data

Surface velocity current fields were obtained from the EuskOOS HF radar station composed by two antennas located at Matxitxako and Cape Higer and covering the SE Bay of Biscay since 2009, a range up to 150 km from the coast. The EuskOOS HF radar is part of JERICO-RI and it is operated following JERICO-S3 project best practices, standards, and recommendations (see Solabarrieta et al., 2016; Rubio et al., 2018, for details). Data consist of hourly current fields with a 5 km spatial resolution obtained from using the gap-filling OMA methodology (Kaplan and Lekien, 2007; Solabarrieta et al., 2021). In total, 85 OMA modes, built setting a minimum spatial scale of 20 km and applied to periods with data from the two antennas, were used to provide maximum spatiotemporal continuity in the HF radar current fields, which is a prerequisite of performing accurate Lagrangian simulations. The application of the OMA methodology has been validated for the Lagrangian assessment of coastal ocean dynamics in the study area by Hernández-Carrasco et al. (2018). HF radar velocities were quality controlled using procedures based on velocity and variance thresholds, signal-to-noise ratios, and radial and total coverage, following standard recommendations (Mantovani et al., 2020). Data subsets were built for the Lagrangian simulations avoiding periods with temporal gaps (still present in the case of the failure of one or two antennas) of more than a few hours. Hourly ERA5-U10-wind fields were obtained from the atmospheric reanalysis computed using the IFS model of the European Center for Medium-Range Weather Forecast (ECMWF) (see C3S, 2019 for details). The ERA5 atmospheric database covers the Earth on a 30 km horizontal grid using 137 vertical levels from the surface up to a height of 80 km and provides estimates of a large number of atmospheric, land, and oceanic climate variables, currently from 1979 to within 3 months of real time. Both HF radar current observations and wind data cover the drifter's emission periods and the selected week-long periods between 2009 and 2021 for riverine litter simulations (see Appendix C for the selected periods).

3.4 Particle transport model

The application of the transport module of the TESEO particle-tracking model (Abascal et al., 2007, 2017a, b; Chiri et al., 2020) was two-fold: (1) to simulate the transport and fate of floating marine litter entering from rivers and reaching the open waters of the SE Bay of Biscay and (2) to estimate a windage coefficient by calibrating the model according to the low-cost buoy trajectories. This module allows for simulating passive particles driven by surface currents, wind, and turbulent diffusion. Particle trajectories were calculated using the following equation:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}\mathit{x}i}{\mathrm{d}t}}={\mathit{u}}_{\mathrm{a}}({\mathit{x}}_i,t)+{\mathit{u}}_{\mathrm{d}}({\mathit{x}}_i,t),\end{array}$$

where u_a and u_d are the advective velocity and diffusive velocity, respectively, for the x_i point and t time. The advective velocity is calculated as the lineal combination of the wind and currents according to

$$\begin{array}{}\text{(2)}& {\mathit{u}}_{\mathrm{a}}={\mathit{u}}_{\mathrm{c}}+C_{\mathrm{d}}{\mathit{u}}_{\mathrm{w}},\end{array}$$

where u_c is the surface current velocity, u_w is the wind velocity at 10 m over the sea surface, and C_d is the wind drag coefficient. The turbulent diffusive velocity is obtained using Monte Carlo sampling in the range of velocities [−u_d, u_d], which are assumed to be proportional to the diffusion coefficients (Hunter et al., 1993; Maier-Reimer and Sündermann, 1982). For each time step Δt, the velocity fluctuation is defined as

$$\begin{array}{}\text{(3)}& \left|{\mathit{u}}_{\mathrm{d}}\right|=\sqrt{{\displaystyle \frac{\mathrm{6}D}{\mathrm{\Delta }t}}},\end{array}$$

where D is the diffusion coefficient, whose value is 1 m² s⁻¹ in accordance with previously modelling work for floating marine litter (Pereiro et al., 2019; Ruiz et al., 2022a). Simulations were forced by HF radar surface current velocity and wind data and interpolated at the particle's position for integrating the trajectories. Beaching along the coast was implemented by a simple approach: if the particle reaches the shoreline, it is identified as beached, and it is removed from the computational process. TESEO has been calibrated and validated by comparing virtual particle trajectories to observed surface drifter trajectories at regional and local scale (Abascal et al., 2009, 2017a, b; Chiri et al., 2019). TESEO is a 3D numerical model conceived to simulate the transport and degradation of hydrocarbons, but it has also been successfully applied to the study of transport and accumulation of marine litter in estuaries (Mazarrasa et al., 2019; Núñez et al., 2019) and in open waters (Ruiz et al., 2022a).

3.4.1 Wind drag estimation

Two simulation strategies were combined for (1) estimating the wind drag coefficient and (2) studying the seasonal behaviour of floating items in the area (Sect. 3.5.2) The wind drag coefficient (C_d) was determined by comparing the observed trajectories provided by the low-cost buoys and the modelled trajectories performed with TESEO. The test was done through different parameterizations of the wind drag coefficient ranging from 0 % to 7 % (Table 2). This range was chosen based on previously floating marine litter studies coupling Lagrangian modelling and observations from satellite drifting buoys (Carson et al., 2013; Stanev et al., 2019; Van Der Mheen et al., 2019). The coefficient providing the lowest error was considered the best coefficient to simulate highly buoyant litter. Due to the grid limitations of the surface currents and wind data in the coastal area, the comparison was not initialized at the launching position of the low-cost buoys (river mouths), but instead it was initialized at the closest grid element that contained valid currents and wind data (Table 1). Observed positions were interpolated into a uniform 1 h time, fitting the metocean temporal resolution. A release of 1000 virtual particles was performed every 4 h at the corresponding observed position (Table 2). Particles were tracked over a 24 h period and the trajectory of the centre of mass of all the particles was computed at every time step to represent the track of the particle cloud. Observations were compared to modelled trajectories using the simple separation distance, which is the difference between the observed and the computed position of the centre of mass at a time step t. The mean separation distance $\stackrel{\mathrm{‾}}{D\left(t^{\mathrm{mod}}\right)}$ was calculated for every modelled position based on the simple separation distance following Eq. (4):

$$\begin{array}{}\text{(4)}& \stackrel{\mathrm{‾}}{D\left(t^{\mathrm{mod}}\right)}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\left|{\mathit{X}}^{\mathrm{mod}}\left(t^{\mathrm{mod}}\right)-{\mathit{X}}^{\mathrm{obs}}\left(t^{\mathrm{obs}}\right)\right|,\end{array}$$

where X^mod (t^mod) and X^obs(t^obs) are the modelled and observed trajectories for the simulation period i of a total of N periods. A mean separation distance curve was computed for every wind drag coefficient derived from the mean separation distance curves of the four buoys. The area beneath the mean separation distance curve was calculated to select the more suitable wind drag coefficient. The area $\stackrel{\mathrm{̃}}{D}$ was calculated as a numerical integration over the forecast period via the trapezoidal method following Eq. (4). This method approximates the integration over an interval by breaking the area down into trapezoids with more easily computable areas:

$$\begin{array}{}\text{(5)}& \stackrel{\mathrm{̃}}{D}\approx \underset{t^{\mathrm{mod}=\mathrm{1}}}{\overset{t^{\mathrm{mod}=\mathrm{24}}}{\int }}\stackrel{\mathrm{‾}}{D\left(t^{\mathrm{mod}}\right)}\mathrm{d}t.\end{array}$$

3.4.2 Lagrangian seasonal simulation of riverine litter items

Seasonal simulations were run for low- and high-buoyancy items to assess the seasonal differences in the transport and fate of floating riverine litter once it has reached the open waters of the SE Bay of Biscay. Particles were released around 2.5 nautical miles off the shoreline due to the complexity in resolving small-scale processes of floating riverine and marine litter behaviour in and close to the river mouths. As parameterizations concerning wind effect linked to the object characteristics are scarce, the optimal wind drag coefficient estimated for the buoys (see Sect. 3.5.1) was accounted for by simulating the behaviour of the objects highly exposed to wind. No wind drag parameterization (C_d = 0 %) was applied for low-buoyancy objects not subjected to the wind effect. A total of 10 periods per season uniformly distributed within the study period (2009–2021) were considered for running the simulations based on the availability of HF radar surface current datasets (Appendix C). In total, 80 simulations (40 for C_d = 0 % and 40 for C_d = 4 %) were run for 7 d. For each simulation, 4000 particles were released in eight rivers (500 per river) assuming that river discharges are equal despite the seasonal variations and the morphological differences between rivers (Table 2). The total number of particles modelled for C_d = 0 % was the same as C_d = 4 %. Post-processing was carried out to compute by river (1) the particles' evolution over the time from their release until their arrival at the shoreline and (2) the particles' distribution on the shoreline, counting the number of beached particles per kilometre of shoreline and indicating the spatial concentration per region.

Table 2 Simulation, release, and physical parameter values for wind drag estimation and floating riverine litter simulations.

4 Results

4.1 Riverine litter characterization

In total 1576 items and 11.597 kg of floating riverine litter were sampled and characterized (Fig. 3). Plastic was the most common type of riverine litter in terms of the number of items (95.1 %) and in weight (67.9 %); they were also frequent glass/ceramics (16.1 %) and cloth/textile items (6.9 %) when counted by weight. The top 10 litter items accounted for 93.3 % by number and 72.6 % by weight of the total riverine litter (Table 3). Plastic/polystyrene pieces between 2.5 and 50 cm and other plastic/polystyrene identifiable items (e.g. food labelling) were the most abundant in terms of number (71.2 %) and weight (16.9 %). Weakly buoyant items encompassed almost 91 % by number and 68 % by weight of litter items (Fig. 4).

Figure 3 Composition of riverine litter by type of material in terms of the number of items and weight. Items were collected by the barrier placed in the Deba River (Gipuzkoa) between April and June 2018.

Table 3 Top 10 (X) riverine litter items collected from the barrier located in the Deba River (Gipuzkoa) between April and June 2018. Items have been ranked by abundance (left) and weight (right) according to the master list categories of beach litter items and classified based on their exposure to the wind effect.

Figure 4 Riverine litter classification based on the exposure to the wind effect. Items were collected from the barrier located in the Deba River (Gipuzkoa) between April and June 2018.

4.2 Wind drag coefficient for drifting buoys

Total distances covered by drifting buoys ranged from 62 to 118 km (Table 1), and they all scattered over the HF radar coverage area. Buoys provided their position data over 385 h before beaching on the Landes and Gipuzkoa shorelines. When compared with numerical trajectories obtained using different C_d parameterizations, the mean separation distance ($\stackrel{\mathrm{‾}}{D\left(t^{\mathrm{mod}}\right)}$ increased nearly linearly with time for all the parameterizations, achieving a maximum separation of almost 14 km at 24 h for C_d = 0 % (Fig. 5). Overall, using no windage parameterization provided the largest $\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{D}}$. Simulations parameterized with C_d = 4 % provided the best results with an average ± standard deviation (SD) of 3.2 ± 1.25 km and a maximum value of 4.85 km at 24 h. When assessing the mean separation distance for all the modelled positions at every observed position of the buoys, the most common range separation distance for C_d = 4 % was 2–4 km (Fig. 6). Hence, a wind drag coefficient of 4 % was applied in the remaining analysis to estimate the behaviour of highly buoyant items.

Figure 5 Mean separation distance between modelled and observed trajectories for each wind drag coefficient. The dark line is the mean curve used for the trapezoidal integration.

Figure 6 Spatial mean distance between modelled and observed trajectories of buoy A, B, C, and D with a drag coefficient C_d = 4 %. Particle trajectories were simulated during 24 h, with a re-initialization period every 4 h. The modelled trajectories are shown in solid lines. Circles represent the mean separation distance at the observed position for all the modelled positions.

4.3 Seasonal trends in floating riverine litter transport and fate

Particle concentration on the shoreline varied between 0 and 258.46 particles km⁻¹ (Fig. 7). Particles parameterized with C_d = 4 % drifted faster towards the coast, notably during the first 24 h. The highest concentrations (>200 particles km⁻¹) were recorded during summer in Pyrénées-Atlantiques for C_d = 4 %, probably due to the seasonal retention patterns within the study area (Appendix D). Although less intensely, C_d = 4 % also led to a high particle concentration in Pyrénées-Atlantiques (106.86 particles km⁻¹) and Gipuzkoa (166.1 particles km⁻¹) during winter. The lowest concentrations (0–20 particles km⁻¹) were recorded for C_d = 0 % after the first 24 h of simulation, particularly during autumn. Overall, Bizkaia was the less impacted region for both windage coefficients (<40 particles km⁻¹). During summer, over 97 % of particles parameterized with C_d = 4 % beached after 1 week of simulation (Fig. 8). In autumn this value fell to 54 %. In contrast, beached particles parameterized with C_d = 0 % were less abundant by the end of the simulations, particularly during spring with less than 25 % of particles trapped in the shoreline.

Overall, the average of particles parameterized with C_d = 0 % was higher when comparing to C_d = 4 % (Fig. 9). Particles released in French rivers and parameterized with C_d = 0 % were less abundant during summer, though this trend was reversed in autumn. For C_d = 0 %, the number of particles released in the Bidasoa River during summer were the least abundant after 1 week of simulation (<200 particles on average). The vast majority of particles released in the Urumea River during winter were floating in the study area by the end of the simulations (479 particles on average). Particles parameterized with C_d = 4 % beached faster during the first 48 h, mainly in summer and for those particles released in the French rivers. During this season, the average number of particles floating in the study area by the end of the simulation ranged between 0 and 250. Similar trends were observed within the same season between rivers, probably influenced by the vicinity of rivers and the spatiotemporal resolution of forcings. Over 40 % of the total particles parameterized with C_d = 4 % and almost 12 % parameterized with C_d = 0 % beached in Gipuzkoa (Fig. 10). During spring, almost 60 % of beached particles parameterized with C_d = 0 % were located Bizkaia. For C_d = 0 %, particles released during summer in the rivers located in the western area of Gipuzkoa drifted longer distances and reached the Landes shoreline. This trend changed during winter, when the vast majority of particles released in Gipuzkoa rivers beached mainly in Gipuzkoa and Bizkaia. Beached particles parameterized with C_d = 0 % experienced more seasonal variations derived from the surface current circulation patterns within the SE Bay of Biscay. For C_d = 4 %, particles beached in Gipuzkoa ranged between 51 % in spring and 38 % in winter, and Bizkaia was the less affected region despite the season. Overall, all regions were highly affected by rivers within or nearby the region itself.

Figure 7 Particle concentration on the Bizkaia, Gipuzkoa, Pyrénées-Atlantiques, and Landes shoreline. The seasonal distribution is shown for C_d = 0 % and C_d = 4 % after 24 and 168 h of simulation.

Figure 8 Seasonal numbers of beached particles parameterized with C_d = 0 % and C_d = 4 % after 168 h of simulation.

Figure 9 Temporal evolution of the particles parameterized with C_d = 0 % and C_d = 4 % throughout the different seasons. The curves represent the average number of particles floating in the water surface by river and for every time step.

Figure 10 Seasonal analysis of the beached particles parameterized with C_d = 0 % and C_d = 4 % per region and river by the end of the simulation period. The nodes of the region correspond to the number of beached particles. The country to which each river belongs – France (Pyrénées-Atlantiques) and Spain (Gipuzkoa) – is shown on the left side of each figure. The width of the node depicts the sum of the beached particles, and the links represent the number of particles beached per river.

5 Discussion

5.1 Riverine litter composition

An artisanal net placed at the mouth of the Deba River enabled sampling riverine litter in the study area during spring 2018. Short and narrow rivers prevail in the SE Bay of Biscay, affected by a strong tidal regime and very intense, stationary and persistent storms (Ocio et al., 2015). Studies aiming at reporting the abundance and composition of floating riverine litter in European rivers date back less than 10 years, and they were performed in larger and more abundant rivers than the Deba River. Despite the morphology and hydrological differences, plastic was the predominant material in the Deba River, as in the Siene (Gasperi et al., 2014), Danube (Lechner et al., 2014) or Rhine River (van der Wal et al., 2015). Similarities were also found when comparing the top 10 list of riverine litter items to rivers located in the north-east Atlantic region. Plastic/polystyrene pieces between 2.5 and 50 cm (71.2 %) top the list in terms of the number of items, and their abundance was slightly higher when compared to north-east Atlantic rivers (54.53 %) (Bruge et al., 2018; González-Fernández et al., 2018). Lower abundances were observed in the Mediterranean (25.01 %) and the Black Sea (13.74 %). Riverine litter items trapped on vegetation or deposited on the riverbank can be degraded by weather conditions (rain, wind, etc.) favouring the fragmentation in plastic pieces before their arrival in the coastal and marine environment (Chamas et al., 2020). The fragmentation can be also influenced by the material and the shape of the litter items (Woods et al., 2021). Differences in plastic/polystyrene pieces between 2.5 and 50 cm abundances can be attributed to a faster fragmentation due to the variations in weather conditions between river basins. However, more detailed analyses on the physical characteristics of litter items (i.e. polymer type) are necessary to fully assess their impact on the occurrence of fragmented plastic pieces. Results are also in line with the ranking list of the top 10 beach litter items across the north-east Atlantic region revealing that single-use plastics (i.e. food containers, bottles, and other packaging) are among the most abundant riverine litter items together with plastic fragments (European Commission, 2018). These results differed from the analysis performed in sea small-scale convergence areas of floating marine litter (“litter windrows”) on the coastal waters of the SE Bay of Biscay, where fishing-related items were the second most abundant sub-category in terms of number after plastic/polystyrene pieces between 2.5 and 50 cm (Ruiz et al., 2020). Substantial differences also exist between riverine litter sampled in the Deba River and floating marine litter assessed by visual observation from research vessels in open waters of the Bay of Biscay (Ruiz et al., 2022a). Differences might be related to the monitoring method and, also, to the size of the items, since small items, such as plastic pieces, can be overlooked by the observer when the visual counting method is applied, contrary to riverine litter samplings for later analysis in the lab. Overall, riverine litter data acquisition is mainly focused on the floating fraction, and the litter loads under the surface water are often ignored. Increasing the quantity of rivers sampled, the frequency, and the riverine water compartments is necessary to establish the composition and trends of riverine litter in the SE Bay of Biscay.

5.2 Wind drag estimation

One of the largest uncertainties for predicting floating riverine and marine litter behaviour is the proper quantification of a wind drag coefficient. Wind drag estimations conducted so far for floating marine litter items range between 0 % and 6 % (Ko et al., 2020; Critchell and Lambrechts, 2016; Neumann et al., 2014) with an upper limit of 10 % (Yoon et al., 2010). However, only a few of them have been validated using observational data (Maximenko et al., 2018; Callies et al., 2017). In this study, data provided by low-cost buoys combined with surface current measurements by HF radar were used as a proxy for modelling the drift of floating litter objects with similar buoy characteristics (density, size, and shape). Results demonstrated that C_d = 4 % was the optimal wind drag coefficient for accurately represent the pathways of the low-cost buoys in the study area. This value can be consistent with the estimations of the partially emerged Physalia physalis for the Bay of Biscay (Ferrer and Pastor, 2017), but it is almost 3 times higher than the maximum wind drag coefficient reported in the area by Pereiro et al. (2018). This can be explained by the fact that buoys used in the experiment remained submerged beneath the sea surface and were less exposed to the wind effect. The estimated wind drag coefficient was also greater than C_d = 3 % observed for the Prestige oil spill accident (Abascal et al., 2009; Marta-Almeida et al., 2013). Indeed, oil spill studies refer to a range of wind drag coefficient between 2.5 % to 4.4 % of the wind speed, with a mean value of 3 %–3.5 % (e.g. ASCE, 1996; Reed et al., 1994). Object characteristics may change over time due to the exposure to wind, waves, UV radiation, seawater, and the attachment of organic material (Kooi et al., 2017; Min et al., 2020). Objects become breakable, and biofouling increases their density, overcoming the positive buoyancy and impacting their trajectory. Investigations so far pinpointed longer timescales (weeks to months and longer) than considered in this study (days) for a significant change on the behaviour of floating objects (Ryan, 2015; Fazey and Ryan, 2016). Consequently, physical variations in the buoy properties were not accounted for the wind drag estimation. The separation distance between observed and modelled trajectories has been commonly used to evaluate the skill of particle-tracking models (Callies et al., 2017; Haza et al., 2019; Aksamit et al., 2020; Abascal et al., 2012). In this study, the purpose was not to evaluate the model accuracy but to estimate the wind drag coefficient for the low-cost buoys. However, the novel approach proposed by Révelard et al. (2021) may be of particular interest for future experiments oriented towards assessing the wind drag coefficient of highly buoyant items drifting for short time periods in the coastal area.

5.3 Seasonal riverine litter distribution by region

It is broadly accepted that the SE Bay of Biscay is polluted with floating marine litter discarded or lost in the marine and coastal area but also with litter originating inland and transported via rivers and run-off. However, detailed studies on riverine litter contribution are still scarce, and modelling efforts combining observations and physical parameterizations of floating litter properties are non-existent. This study shows that the exposure to the wind effect largely controls the transport and coastal accumulation of floating marine litter in the SE Bay of Biscay, with concentrations varying between regions and over time. Concentrations in Pyrénées-Atlantiques and Gipuzkoa differed widely from the other studied regions. Indeed, the highest concentrations occurred in both regions during summer for low- (100–120 particle km⁻¹) and high-buoyancy items (> 200 particles km⁻¹). A higher number of particles beached in Gipuzkoa during summer when compared to Pyrénées-Atlantiques, but concentrations were lower since the Basque shoreline is longer. The pathways and fate of low-buoyancy items reflect the seasonal surface water circulation patterns in the SE Bay of Biscay. Results are in line with findings provided by Declerck et al. (2019) who pinpointed a higher coastal retention in the area during spring and summer. Weakly buoyant objects remained floating at the coastal waters and highly buoyant objects tended to beach remarkably faster as reported in literature by Rodríguez-Díaz et al. (2020). However, long-term data collected by in situ observations of beached litter across the different regions are necessary to validate the large seasonal variations and to assess the reliability of concentration levels for addressing riverine litter issue in priority regions with heavily polluted coastlines.

5.4 Rivers as key vectors of riverine litter

The interpretation of the spatial and temporal riverine litter distribution by river can be challenging since riverine litter fluxes in the study area are highly uncertain. In the study area, two major assumptions were made regarding the river systems: (1) the same river discharge for all rivers and (2) the same river discharge for all seasons. This means that the same amounts of riverine litter were allocated for every river regardless of the differences in the width and depth and the seasonal flow variations. Since each river basin has its own particularities, future modelling approaches should be adapted to the morphology and hydrological conditions of the catchment area. Other drivers, such as the land use or population density, can be a determining factor for the amount of mismanaged litter that could contribute to riverine litter fluxes (Schmidt et al., 2017; Schuyler et al., 2021). It is also necessary to further investigate if higher river flows in the area are directly related to an increased discharge of riverine litter since analysis already performed in different river basins shows contradicting relations between the occurrence of riverine litter and river fluxes (van Emmerik and Schwarz, 2020). Along with the complex nature of qualifying riverine litter fluxes, litter behaviour in the coastal area of the SE Bay of Biscay is still in its early stage, and much has yet to be revealed. Particular attention should be paid to Pyrénées-Atlantiques and Gipuzkoa, as the main impacted regions in the studied area. Rivers in the study area are mainly located in Gipuzkoa, which favours the accumulation of floating litter in this region regardless of the season. Regional coordination should be reinforced due to the transboundary movement of floating riverine litter in the study area and reasonable efforts oriented towards retaining or removing riverine litter as clean-up measures in the riverbanks should be investigated to avoid litter being transported to the coastal and marine environment.

5.5 Model limitations

The interaction between floating litter and the shoreline is highly complex and relies on many processes including waves and tides. Indeed, waves and tides can constrain coastal accumulation since they can resuspend and transport litter back into the ocean (Brennan et al., 2018; Compa et al., 2022). The geomorphology can also affect the retention of litter washing ashore. Sandy beaches tend to be more efficient at trapping and accumulating litter than rocky areas, which favour litter fragmentation (Robbe et al., 2021; Weideman et al., 2020). How these processes contribute to the actual beaching is unknown, and they cannot be resolved yet at a suitable resolution (Melvin et al., 2021). In this study, particles were released in open waters, and once they reached the shoreline, they were classified as beached. The tidal effect and the wave-induced Stokes drift were not accounted for to avoid introducing more uncertainties. However, further field and laboratory experiments to better understand how these processes influence floating litter behaviour on the coastline are recommend. For future research, it is also important to consider exploring the effect of the type of shoreline on coastal accumulation. In this study, a constant diffusion coefficient of 1 m² s⁻¹ was regarded as a pragmatic choice based on previous modelling work for floating marine litter. However, more field measurements are necessary to accurately assess the influence of the diffusion process on the transport of floating marine litter.

6 Conclusions

The SE Bay of Biscay has been described by global and regional models as an accumulation zone for floating marine litter. However, detailed studies on floating riverine litter behaviour once items arrive in open waters are still scarce. Based on HF radar current observations and a wind dataset, this contribution tries to fill this gap by providing insights into how low- and high-buoyancy litter released by several rivers of the SE Bay of Biscay may affect the nearby regions seasonally in terms of concentration and beaching. Analysis of riverine litter samples collected by a barrier placed in the study area showed that low-buoyancy objects were predominant, although highly buoyant objects were also relevant in terms of weight. Simulations for assessing the seasonal trends of floating riverine litter transport and fate were performed with the Lagrangian model TESEO. To properly integrate the differences in litter buoyancy, simulations were parameterized with a wind drag coefficient for low- and high-buoyancy items. The wind drag for highly buoyant items was estimated by comparing the observed and the modelled positions of four drifters. The developed low-cost buoys proved to be suitable to provide real-time trajectories of highly buoyant objects exposed to wind. However, drifters with different characteristics should be used in future studies to account for the windage effect on different types of items. The transport and fate of both high- and low-buoyancy items released by rivers was calculated by season. Highly buoyant items rapidly beached (in less than 48 h), particularly in summer and winter; in contrast, despite the season over two-thirds of low-buoyancy items remained floating after 1 week of being released. This highlights the discrepancy between the behaviour for low- and high-buoyancy objects and the importance of parameterizing the windage effect in order to accurately predict riverine litter accumulation in the coastal area of the SE Bay of Biscay. Beached particles were mainly found in Gipuzkoa regardless of the season and the wind drag coefficient. Overall, the less affected region was Bizkaia with the exception of a spring period for low-buoyancy items. Despite the season, most of the riverine litter remained in the study area and rivers polluted the regions within the river basin or surrounding it. Investigating what beaches are most likely to accumulate large quantities and the contribution per river can provide relevant input to response operations after storm events in the short to medium term and can also support the identification of priority rivers for a monitoring programme, assisting adapted intervention of riverine pollution regionally in the future.

Appendix A: Floating barrier for riverine litter collection

Figure A1 Floating barrier (a) and installation in the Deba River (Gipuzkoa) (b).

Appendix B: Riverine litter classification based on the exposure to the wind effect

Table B1 Data were gathered from surveys carried out during spring 2018 in the Deba River (Gipuzkoa).

Appendix C: Selected seasonal week-long periods from the HF radar (2009–2021)

Table C1 Periods selected between 2009 and 2021 based on the availability surface current datasets provided by the HF radar.

Appendix D: Seasonal mean current and wind fields (2009–2021)

Figure D1 Mean current (a) and wind fields (b) in the study area during each season for the selected periods between 2009 and 2021. The colour bars represent the magnitude of current and wind speed. The arrows indicate the current and wind mean direction and are scaled with currents and wind speed (data sources: HF radar – EuskOOS, https://www.euskoos.eus/en/data/basque-ocean-meteorological-network/high-frequency-coastal-radars/; ERA5, https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5, last access: 11 November 2019).