Ocean Plastic Assimilator v0.2: assimilation of plastic concentration data into Lagrangian dispersion models

A numerical scheme to perform data assimilation of concentration measurements in Lagrangian models is presented, along with its first implementation called Ocean Plastic Assimilator, which aims to improve predictions of the distributions of plastics over the oceans. This scheme uses an ensemble method over a set of particle dispersion simulations. At each step, concentration observations are assimilated across the ensemble members by switching back and forth between Eulerian and Lagrangian representations. We design two experiments to assess the scheme efficacy and efficiency when assimilating simulated data in a simple double-gyre model. Analysis convergence is observed with higher accuracy when lowering observation variance or using a circulation model closer to the real circulation. Results show that the distribution of the mass of plastics in an area can effectively be improved with this simple assimilation scheme. Direct application to a real ocean dispersion model of the Great Pacific Garbage Patch is presented with simulated observations, which gives similarly encouraging results. Thus, this method is considered a suitable candidate for creating a tool to assimilate plastic concentration observations in real-world applications to estimate and forecast plastic distributions in the oceans. Finally, several improvements that could further enhance the method efficiency are identified.

1 Introduction

Plastic pollution reveals itself to be an urgent matter if humans are to preserve their oceans. Previous publications such as Lebreton et al. (2018) reviewed how plastics are rapidly accumulating in the oceans and concentrating in oceanic gyres. As public and private ventures set out cleanup goals, accurate and regular forecasts of the state of plastics in the oceans are becoming necessary.

A modeling framework is currently undergoing development at The Ocean Cleanup towards this goal, as the company set itself out to clean 90 % of the oceans' floating macroplastics by 2040. It is used to assess and improve our ability to perform the largest cleanup in history.

This framework, the results of which are presented in Lebreton et al. (2018), is built upon the Pol3DD Lagrangian dispersion model and presented in Lebreton et al. (2012). In this model, virtual particles representing plastics are generated and let drift over time using current data extracted from the oceanic circulation modeling system HYCOM (HYbrid Coordinate Ocean Model; see Bleck, 2002). Results from this model are compared with two other plastic forecast models in van Sebille et al. (2015).

While the Lebreton et al. (2012) modeling framework has already produced valuable results, it is not able to assimilate observations and update forecasts accordingly yet. However, as the company prepares to release a number of systems to clean the ocean, it will soon have numerous sources of data-collecting devices measuring the concentration of plastics in the oceans. Therefore, we believe it is timely to develop a method to assimilate incoming real-time observations.

Methods to assimilate plastic concentration observations over a Lagrangian dispersion model are in the early development stage (Lermusiaux et al., 2019). However, earlier studies dealing with data assimilation applied to the atmospheric dispersion of particles around polluting facilities, such as Zheng et al. (2007), have been published.

This paper introduces Ocean Plastic Assimilator v0.2, a numerical scheme developed to assimilate plastic concentration data into 2D Lagrangian dispersion models. Section 2 formulates the method, and Sect. 3 then describes its initial implementation and application. For this proof-of-concept paper, we use a dispersion simulation generated with the OpenDrift framework in a controlled environment based on a double-gyre analytical flow field. The assimilation results are presented in Sect. 4. Real-world application perspectives and future developments that could further improve the method are discussed in Sect. 5. Finally, in Sect. 6, we present a direct application of the method to a dispersion model of the actual Great Pacific Garbage Patch, with simulated observations sampled from another simulation.

2 Method

This section formulates our methodology to perform data assimilation of plastic concentration (or density) observations in any 2D Lagrangian dispersion model using an ensemble Kalman filter (EnKF). It includes the two representations of data (Eulerian and Lagrangian) being used for this process, the transformation between Eulerian and Lagrangian space, the ensemble assimilation method itself, and model ensemble initialization.

2.1 Representations of data

The distribution of the mass of plastics in a Lagrangian dispersion model is represented by weighted particles drifting according to a flow field in a 2D domain. Each virtual particle represents a drifting concentration of plastics. In turn, virtual concentration measurements are collected at fixed locations (grid points) within the studied 2D domain, i.e., an Eulerian representation of the plastic mass distribution.

Our method aims to assimilate concentration observations collected in the Eulerian representation and update the Lagrangian representation accordingly. One cycle of this process consists of projecting particle weights on the concentration grid, assimilating observation data into the concentration grid, projecting grid cell concentration updates on particle weights, and finally letting particles drift until the next assimilation time step. This procedure is summarized in Fig. 1.

Figure 1 Schematic depiction of the four steps of our method.

The complete workflow requires the following:

• an assimilation method,

• a dispersion model along with the flow field used,

• projection methods to go back and forth between Eulerian and Lagrangian representations, and

• prior estimates for model parameters and uncertainties.

2.2 Procedure

This section presents our procedure for a set of N_p particles drifting in a gridded domain, with a grid size (m,n) and indices i and j to designate a grid cell. An ensemble Kalman filter works by running different simulations, or ensemble members, simultaneously with variations in model parameters (e.g., initial conditions). We use N_e members in the following.

2.2.1 Projecting weights on densities

At each step t, we define the following:

• ${\mathit{w}}_t^{\mathrm{f}}$ the forecast weights vector of size N_p (kg);

• ${\mathit{x}}_t^{\mathrm{f}}$ the forecast density vector computed after projecting ${\mathit{w}}_t^{\mathrm{f}}$ on the density grid (kg m⁻²);

• y_t the density observation vector (kg m⁻²), with its error covariance matrix R;

• ${\mathit{x}}_t^{\mathrm{a}}$ the analyzed density vector computed by assimilating observations y_t in ${\mathit{x}}_t^{\mathrm{f}}$ via the ensemble Kalman filter (kg m⁻²);

• ${\mathit{w}}_t^{\mathrm{a}}$ the analysis weight vector computed by projecting on ${\mathit{w}}_t^{\mathrm{f}}$ the corrections computed on ${\mathit{x}}_t^{\mathrm{a}}$ (in kg); and

• ${\mathrm{\Delta }}_{i,j,t}$ the set of particles present at step t in grid cell i,j.

To start, for grid cell i,j with area A_i,j, ${\mathit{x}}_t^{\mathrm{f}}$ is computed with the formula

$$\begin{array}{}\text{(1)}& ({\mathit{x}}_t^{\mathrm{f}}{)}_{i,j}={\displaystyle \frac{{\sum }_{p\in {\mathrm{\Delta }}_{i,j,t}}({\mathit{w}}_t^{\mathrm{f}}{)}_p}{A_{i,j}}}.\end{array}$$

In the following, we omit sub-index t when all the operations are performed at the same time step t.

2.2.2 Assimilating with the ensemble Kalman filter (EnKF)

Our assimilation step relies on the use of ensemble Kalman filtering, as described in Evensen (2003). This method is derived from Kalman filtering and is notably suitable for situations in which the model is not an easily invertible matrix (used in standard Kalman filtering) and one cannot efficiently compute an adjoint (used in extended Kalman filtering).

Standard Kalman filtering allows computing the analysis state using a single equation. In standard Kalman filtering, the forecast state vector x^f (in this case, the densities) and the analysis vector x^a are linked with

$$\begin{array}{}\text{(2)}& {\mathit{x}}^{\mathrm{a}}={\mathit{x}}^{\mathrm{f}}+\mathbf{K}(\mathit{y}-\mathbf{H}{\mathit{x}}^{\mathrm{f}}).\end{array}$$

H is the observation matrix that maps the state x^f to the observation space of y.

The Kalman gain matrix K is defined by the following equation.

$$\begin{array}{}\text{(3)}& \mathbf{K}={\mathbf{P}}^{\mathrm{f}}{\mathbf{H}}^T({\mathbf{HP}}^{\mathrm{f}}{\mathbf{H}}^T+\mathbf{R}{)}^{-\mathrm{1}}\end{array}$$

R is the observation error covariance matrix. P^f is the forecast error covariance matrix. When using a Kalman filter, P^f is in principle meant to be computed from the previous state by application of the forward integration matrix operator, but this is generally too computationally expensive and impractical. Here, we use ensemble Kalman filtering, whereby the P^f matrix computation is approximated by relying on an ensemble of simulations.

Ensemble members are different instances of our simulation with different initializations. For ensemble member $k\in \left[\right|\mathrm{1},N_{\mathrm{e}}\left|\right]$, we write ${\mathit{x}}_k^{\mathrm{f}}$ the forecast state vector and ${\stackrel{\mathrm{‾}}{\mathit{x}}}^{\mathrm{f}}$ the ensemble average:

$$\begin{array}{}\text{(4)}& {\stackrel{\mathrm{‾}}{\mathit{x}}}^{\mathrm{f}}={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{e}}}}\sum _{k=\mathrm{1}}^{N_{\mathrm{e}}}{\mathit{x}}_k^{\mathrm{f}}.\end{array}$$

Accordingly, the computation of P^f can be accomplished using the formula

$$\begin{array}{}\text{(5)}& {\mathbf{P}}^{\mathrm{f}}={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{e}}-\mathrm{1}}}\sum _{k=\mathrm{1}}^{N_{\mathrm{e}}}({\mathit{x}}_k^{\mathrm{f}}-{\stackrel{\mathrm{‾}}{\mathit{x}}}^{\mathrm{f}})({\mathit{x}}_k^{\mathrm{f}}-{\stackrel{\mathrm{‾}}{\mathit{x}}}^{\mathrm{f}}{)}^T.\end{array}$$

Each ensemble member k is then updated using Eq. (2) with x_k instead of x.

2.2.3 Projecting the density updates on particles

Several ways of projecting the density updates (step 3 in Fig. 1) can be thought of. In the Ocean Plastic Assimilator v0.2, we simply choose to update the weights by uniformly distributing the density correction ratio of a grid cell i,j among the particles in the same box using this formula:

$$\begin{array}{}\text{(6)}& \forall p\in {\mathrm{\Delta }}_{i,j},({\mathit{w}}^{\mathrm{a}}{)}_p={\displaystyle \frac{({\mathit{x}}^{\mathrm{a}}{)}_{i,j}}{({\mathit{x}}^{\mathrm{f}}{)}_{i,j}}}({\mathit{w}}^{\mathrm{f}}{)}_p.\end{array}$$

In this equation, (x^f)_i,j cannot be null when a grid cell i,j contains particles (see Eq. 1) unless all particles have null weights. While extremely unlikely (we did not encounter this phenomenon during our numerous tests), particles with exactly null weights have to be taken out of the simulation.

This heuristic was chosen primarily for its simplicity and its computational efficiency. The multiplicative approach also tends to prevent computing negative weights if the density analysis is lower than the density forecast.

Finally, for step 4 in Fig. 1, since the dispersion model changes particles' positions but not their weights when integrating, the forecast weights at time t+1 are

$$\begin{array}{}\text{(7)}& {\mathit{w}}_{t+\mathrm{1}}^{\mathrm{f}}={\mathit{w}}_t^{\mathrm{a}}.\end{array}$$

2.2.4 Initialization

As stated by Evensen (2003) the ensemble Kalman filter requires the initial ensemble to sample the uncertainty in variables that we want to update with data assimilation. In this article, we focus on our method's ability to compute the correct total mass of particles drifting. For this reason, we normally distribute the members' initial total masses with a mean μ_e and standard deviation σ_e. If we write M_k the initial total mass for ensemble member k, we thus have

$$\begin{array}{}\text{(8)}& M_k\sim N({\mathit{\mu }}_{\mathrm{e}},{\mathit{\sigma }}_{\mathrm{e}}).\end{array}$$

Finally, we attribute an initial weight of $M_k/N_{\mathrm{p}}$ to each particle.

3 Implementation and test-case setup

This section presents the Python implementation of the aforementioned method, called Ocean Plastic Assimilator (v0.2). We then describe the Lagrangian dispersion model (OceanDrift) used to generate double-gyre dispersion simulations and the experiments created with it to observe how our method performs in a controlled environment.

3.1 Python implementation of the Ocean Plastic Assimilator

This first implementation is coded in Python (see Peytavin, 2021a, for the repository). It is meant as a stand-alone program using dispersal model output data as input, formatted as a NetCDF4 dataset containing particle coordinates in a given space and time domain, along with their weights. It is assumed that the advection in the dispersion model does not depend on particle masses. In the more general case, one would have to run the model again after each assimilation time step, as a change of a particle mass could change its future trajectory.

Once loaded, the input weights are duplicated in N_e arrays, and the program runs the assimilation scheme presented in the previous section in a time loop, taking observations from an input data frame at each time step. The assimilator can also take one additional dispersion simulation output from which it samples observations to assimilate at each time step. This is the approach used in the following test case.

This implementation leverages the use of arrays and the fact that we only use one simulation for all ensemble members to perform vectorized computations for the computation of P^f in Eqs. (1) and (6). It also allows computing ${\mathrm{\Delta }}_{i,j,t}$ only once for all ensemble members. Some parts of the algorithm are also executed with the just-in-time compiler Numba (see Lam et al., 2015) in order to run faster.

This implementation allows our algorithm to perform each following test case with repeated assimilation of two observation points during 2000 time steps in a (60,40) gridded domain in less than 30 min on a modern laptop using about 3 GB of storage and 2 GB of RAM.

Running the assimilator on a dispersion output and not inside a dispersion model allows it to work on outputs from different models, as long as the data are appropriately formatted. Future implementations could also offer the option of running online (i.e., embedded inside a dispersion model), which could allow more flexibility and possibilities, as discussed in Sect. 6.2.3.

3.2 Double-gyre plastic dispersion using the OceanDrift model

In order to create our test cases, we first need a dispersion model and a flow field. We chose the OceanDrift model from the Norwegian Lagrangian trajectory modeling framework OpenDrift (see Dagestad et al., 2018). It was chosen mainly for its simplicity and the fact that OpenDrift embeds a module to generate a dispersion based on a 2D double-gyre flow field.

Figure 2 Generated particles (a) in a double-gyre flow field with OpenDrift and the corresponding plastic concentration field in particles per grid area (b). The domain grid size is 60×40.

This field consists of two gyres periodically moving closer then farther away in an enclosed area. It is a simple field but complex enough to stir and disseminate particles and is regularly used as a standard case to study time-varying flows, for example in Guo et al. (2018). The evolving currents are generated using an analytical field¹. The equations generating this 2D, time-varying, deterministic field are as follows.

$$\begin{array}{}\text{(9)}& {\displaystyle }\begin{array}{rl}u& =-{\displaystyle \frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}y}}=-\mathit{\pi }A\sin \left(\mathit{\pi }f\right(x,t\left)\right)\cos \left(\mathit{\pi }y\right)\\ v& ={\displaystyle \frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}x}}=\mathit{\pi }A\cos \left(\mathit{\pi }f\right(x,t\left)\right)\sin \left(\mathit{\pi }y\right){\displaystyle \frac{\mathrm{d}f}{\mathrm{d}x}}\end{array}\text{(10)}& {\displaystyle }\left\{\begin{array}{l}f(x,t)=a\left(t\right)x^{\mathrm{2}}+b\left(t\right)x\\ a\left(t\right)=\mathit{ϵ}\sin \left(\mathit{\omega }t\right)\\ b\left(t\right)=\mathrm{1}-\mathrm{2}\mathit{ϵ}\sin \left(\mathit{\omega }t\right)\end{array}\right.\end{array}$$

The dimensionless domain size for these equations is $[\mathrm{0},\mathrm{2}]\times [\mathrm{0},\mathrm{1}]$.

Parameter A is the circulation amplitude, ω is the frequency of oscillation of the gyres, and ϵ is the amplitude of the gyre oscillation relative to the steady state.

Particles are then generated and advected using the OceanDrift Lagrangian model from the Norwegian trajectory modeling framework OpenDrift (Dagestad et al., 2018). Figure 2 shows such a dispersion and the associated concentration field.

Thus, we can generate different dispersion simulations by changing the initial particle position seed, which changes the distribution of particle trajectories and the initial masses of the particles. We can also change the flow field parameters A, ω, and ϵ.

In the following section, we modify the flow field parameters and the particle position seeds to create assimilation test cases that use two simulations: a reference and a forecast. We then sample observations from the reference simulation and assimilate them inside the forecast simulation. By doing so, we mimic assimilating real concentration data into an uncertain flow field in the presence of model error.

3.3 Assimilation experiment setup

In order to assess and quantify the efficacy of the assimilator in different cases, we designed two experiments.

The first one aims to verify that, when the forecast flow field reproduces the reference flow field accurately, our implemented scheme can correct an incorrect total mass guess. It also intends to check that the estimate gets better when the observation error gets lower, as one would generally expect.

The second experiment aims to assess the assimilator's behavior and efficacy when the forecast flow field is slightly different from the reference by changing the double-gyre parameters A and ϵ.

In both experiments, we run several test cases to assimilate observations taken from a reference simulation into a forecast simulation using the assimilator. Then, we compute the total plastic mass estimation error and the concentration field root mean square error (RMSE) to assess how close the assimilated forecast gets to the reference situation. This procedure is depicted in Fig. 3.

Figure 3 Schematic depiction of a test case using a reference and a forecast simulation.

In each test case, the Ocean Plastic Assimilator is executed over the course of 2000 time steps. The double-gyre size, which is $[\mathrm{0},\mathrm{2}]\times [\mathrm{0},\mathrm{1}]$, is dimensionless; this means that the time step is dimensionless too. However, if the flow field was the size of the Great Pacific Garbage Patch, then with A=0.1 the time step would be of the order of a day.

Over the double gyre, we define a gridded domain of size (60,40) and select two fixed observation points to run each assimilation test case. This sampling pattern can be thought of as representing a set of moorings that one may deploy in the real ocean. H is defined as a matrix that subsets (x)_i,j to two points of observations.

For the ith point, the measurement is simulated by adding a random error to x_i such as

$$\begin{array}{}\text{(11)}& y_i=\text{max}(x_i+N(\mathrm{0},{\mathit{\sigma }}_{\mathrm{rel}}{x}_i),\mathrm{0}).\end{array}$$

To compute matrix R, we choose to model the observation error as a sum of an additive error σ₀ and a multiplicative relative error σ_rel. As such, with y_i the value measured at the ith observation point is

$$\begin{array}{}\text{(12)}& \mathbf{R}=\text{diag}({\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}+({\mathit{\sigma }}_{\mathrm{rel}}{y}_{\mathrm{1}}{)}^{\mathrm{2}},{\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}+\left({\mathit{\sigma }}_{\mathrm{rel}}{y}_{\mathrm{2}}{)}^{\mathrm{2}}\right).\end{array}$$

In the following, unless specified otherwise, we use N_e=10, σ_e=0.05, N_p=25 000, σ₀=0.1, and σ_rel=1 %. The coordinates of the two observation points are the following pairs: $(\mathrm{12},\mathrm{4}),(\mathrm{55},\mathrm{27})$.

4 Results

4.1 Estimating the total plastic mass in the forecast

In this first experiment, we want to assess the ability of our newly implemented scheme to estimate the total mass of plastics in the reference simulation correctly.

First, we generate a reference situation using ϵ=0.25, A=0.1, and $\mathit{\omega }=\mathrm{2}\mathit{\pi }/\mathrm{10}$. We input the same parameters to integrate the particle trajectories in the forecast simulation. By doing so, we are in a position where we understand the flow of the reference situation correctly, but we do not know the total mass of plastics drifting. In the following, M_ref=25 000 is the constant total mass of the reference situation.

We initiate five different forecasts with μ_e=0.25M_ref, 0.5M_ref, M_ref, 2M_ref, and 5M_ref. Observations are collected (and later assimilated) at each time step on two observation points that could, for example, represent a pair of moored instruments.

Figure 4 Evolution of total mass over time for five different forecast simulations with five different initial total masses (Table 1) over 100 assimilation iterations (a) and 2000 iterations (b). The total mass evolution of the reference simulation is indicated by a solid line.

Figure 4 shows the evolution of the forecast total mass for each simulation. Forecasts starting with an initial total mass lower than approximately 0.82M_ref have their total mass rise, while those starting with higher total mass have their total mass fall. Final total plastic masses in the forecast after 1900 steps of assimilation for each simulation are presented in Table 1. Overall, the forecast total masses seem to converge towards a similar value of approximately 0.82M_ref, from which we can conclude that in this situation the method makes an 18 % error.

Table 1 Final total mass (FTM) relative to M_ref and the concentration field RMSE for five different forecast simulations with five different initial total masses μ_e. RMSE_f and RMSE_∅ are the concentration field RMSE at the end of simulations with and without assimilation of observations.

Another indicator of the correctness of a simulation can be computed from the concentration field at each step. For one of the forecasts (with μ_e=2), we analyze the distribution of concentration errors over the gridded domain and through time (Fig. 5). We observe a decrease in the mean absolute percentage error and a decrease in the absolute percentage error standard deviation. We also observe that this distribution does not contain overly large values.

We also compute the concentration field root mean square error (RMSE_f) at the end of the simulation after assimilating and RMSE_∅ at the end of a simulation with no assimilation. Values in Table 1 indicate a clear improvement of the RMSE when the initial total mass was erroneous and a stable one compared to no assimilation when the initial total mass was correct.

Overall, this points to an improvement in the forecast concentration field over time, thanks to data assimilation.

Figure 5 Evolution of the error field between the reference concentration field and the forecast concentration field, in percent, for μ_e=2. At each time step, the error field is computed and the distribution of the absolute errors in each cell, in percent of the cell reference concentration, is depicted in the box plots. Dots outside whiskers represent outliers, and the triangle is the mean.

Finally, in order to assess the method accuracy depending on observation errors, we set μ_e=2 and run simulations with different values of σ_rel. FTM and RMSE are then computed and presented in Table 2.

Table 2 Parameters and metrics for assimilation simulations with different values of σ_rel, with μ_e=2. FTM is the final total mass, and RMSE_f and RMSE_∅ are the concentration field RMSE at the end of simulations with and without assimilating.

We find that decreasing σ_rel increases the final total mass of the forecast, getting it closer to 1, while the RMSE decreases. This demonstrates that the forecast bias can be reduced by decreasing the observation error, as one would usually expect of a data assimilation method.

4.2 Impact of physical model errors

In this second experiment, we change the parameters used to generate the currents of the reference simulation double gyre. For example, the impact of a modification of ϵ on the generated flow field is illustrated in Fig. 6. By assimilating observations from reference situations with different double-gyre parameters, we can observe the effects of having an erroneous physical dispersion model when assimilating data.

We initiate the forecast with an erroneous initial total mass of 2M_ref and expect that the best total mass predictions will arise from assimilation simulations with the closest flow field.

Figure 6 Flow fields at t=2.5 s for two double-gyre simulations with (a) ϵ=0.1 and (b) ϵ=0.5.

The forecast simulation is generated using ϵ_ref=0.25, A_ref=0.1, and ${\mathit{\omega }}_{\mathrm{ref}}=\mathrm{2}\mathit{\pi }/\mathrm{10}$.

We then generate different reference simulations with different values of A and ϵ, and we try assimilating observations sampled from each of them into the forecast.

We find that data assimilation remains effective and that simulations run with values of ϵ and A closer to ϵ_ref and A_ref lead to better estimations of the total mass and concentration field after some time as one might expect (Fig. 7 and Table 3).

Table 3 Parameters and metrics for simulations with different values of A and ϵ for the reference simulation. FTM is the final total mass, and RMSE_f and RMSE_∅ are the concentration field RMSE at the end of simulations with and without assimilating.

Figure 7 Evolution of the total plastic mass in the forecast simulation for five different runs with varying values of double-gyre parameters A and ϵ, along with the total plastic mass in the reference simulation.

This result illustrates that the assimilation method can be robust to unknown model errors.

5 Application to the Great Pacific Garbage Patch

In this section, we present an application to real-world global dispersion models. As before, we sample observations from one simulation and assimilate them into another in order to mimic the assimilation of observations that could be collected daily by a pair of moorings deployed in the real ocean. We just use an estimate of real ocean currents in place of the simplified double gyre defined in Eq. (9).

We generate two global dispersion simulations with the Lagrangian dispersion model presented by Lebreton et al. (2012). In both cases, the circulation model uses output from the HYbrid Coordinate Ocean Model (see Bleck, 2002), available every 6 h at 0.08^∘. This estimate includes Ekman transport and convergence, as well as mesoscale eddies. The first simulation has particles seeded along the coasts of 192 countries depending on reported garbage input estimates. The second simulation has particles seeded at river mouths only based on estimates of their outflow of plastics. Both generation models are described in the Supplement of Lebreton et al. (2018). A model spin-up was done from 1993 to the end of 2011.

We initialize plastic particle masses generated in the coastal-seeded model depending on their release year. If x is the time spent (in fraction of years) since the beginning of the simulation, then $w_{\mathrm{p}}=\mathrm{1}+x+\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\sin \left(\mathit{\pi }\right(\mathrm{2}x+\mathrm{1}\left)\right)$ is the mass of particles, in tonnes, seeded at time x. This formula increments particle masses by 1 t each new release year, with some periodic variability. The particle masses in the river-seeded simulation are initialized to 1 t regardless of their release date. By doing so, we mimic a situation in which we underestimate the yearly increase in plastic mass input to the ocean.

The gridded domain has a resolution of 0.5^∘, with 80 by 44 points, going from 165 to 125^∘ W and from 23 to 45^∘ N. Throughout 2012, we sample two observations per day at positions 152.5^∘ W, 29^∘ N and 140^∘ W, 35^∘ N from the coastal-seeded dispersion simulation and assimilate them in the river-seeded dispersion model. We use N_e=10, σ_e=50, N_p=25 000, σ₀=0.1, and σ_rel=1 %.

Our method is able to predict the total mass of floating plastics with a 17 % error and to divide the concentration field RMSE by 4 (Fig. 8). The computations take about an hour to run on a standard laptop.

The updates to the concentration field are presented in Fig. 9, which shows that, as expected, the assimilated forecast has increased concentrations.

Figure 8 (a) Evolution of total plastic mass in the domain through 2012 for the reference simulation and the forecast simulation. (b) Evolution of the concentration field RMSE in the assimilation domain through the year 2012.

Figure 9 Concentration field updates at the end of the assimilation cycle, with the two observation locations in blue. This field is the difference between the forecast concentration field at the end of the year 2012 with assimilation and the same without assimilation.

Further experimentation will be required to assess the benefits of using this method in real-world use cases with real data. However, these results confirm the potential skill of our method, even in the presence of sizable model error.

6 Discussion and perspectives

6.1 Towards an application to real-world data

In this proof-of-concept paper, we placed ourselves in a controlled environment to assess the efficacy of the method. In the future, our goal will be to eventually apply the method to real data by replacing the simulated reference situation observations with real-world observations, and the previous results can help in understanding what might happen in assimilating real-world data. The fact that replacing the analytic circulation field by a real-world one (in Sect. 5) did not prevent the method from improving the forecast is viewed as an encouraging first step in that direction.

In Fig. 7a and b we observed that the more accurate the underlying dispersion model is, the more accurate the assimilation result is. For our method to be applied successfully to a real global plastic assessment model, its dispersion prediction would have to be accurate enough. Ongoing work, which is focused on identifying model error sources and refining statistical priors, should benefit the planned application to real data (e.g., Maximenko et al., 2012; van Sebille et al., 2020; Meijer et al., 2021).

Conveniently, we observed that the forecast total mass gets higher when the dispersion model is more accurate, thus acting, in a way, like a score. As a result, we might discriminate between dispersion models based on this method's output by selecting the ones that output the highest total mass.

6.2 Future developments

Amongst the potential applications of the presented method, one might highlight the evaluation and design of real observational strategies. Here we considered one hypothetical, albeit plausible, scenario which might represent the deployment of a few relatively accurate moorings. In future studies it would be interesting to investigate how data coverage in space and time may affect forecast skill in more detail, for example, or use this data assimilation system as a benchmark for proposed field campaigns. Several directions to further develop the method and make it more accurate also seem worth considering, as outlined below.

6.2.1 Improving the filter

Throughout the last 2 decades, the ensemble Kalman filter has been extensively developed and improved, with numerous variants published in the scientific literature. Using different ensemble sampling strategies or a square root algorithm was described as a way to improve accuracy in Evensen (2004). Other solutions include inflating the ensemble before assimilating (see Anderson, 2007), resampling the ensemble, or using a method to assimilate observations locally by adding a Schur product with a so-called correlation matrix in the computation of the Kalman gain in Eq. (3) (see Houtekamer and Mitchell, 2002). Assimilating locally around observation locations could also have the advantage of further improving the geography of the concentration field, which would translate in reduced values of RMSE_f.

6.2.2 Decoupling the positions of the particles for all ensemble members

The method presented here uses the same dispersion simulation as a base for the trajectories of the particles for all ensemble members. In all members, the particle positions through time are the same; the only variables that differ are the particle masses. In particular, the particle trajectories are the same in each member. This approach greatly reduces the storage cost and increases computation speed.

However, it significantly lowers the diversity of the ensemble, so in future work one might want to decouple the ensemble member trajectories, i.e., have a unique set of trajectories for each member. We anticipate that extending the method to use an ensemble with diverse particle simulations should help the forecast converge towards a concentration field closer to the reference one. We regard this possibility as a leading candidate to make the method even more accurate.

6.2.3 Studying other projection operators

In Sect. 2.2.3, we presented a simple way to update particle weights after assimilating density observations through Eq. (6). Different possibilities for performing this step have been thought of, some of which we think may be worth investigating further. Another simple approach would be to apply an additive correction instead of the multiplicative correction used in Eq. (6):

$$\begin{array}{}\text{(13)}& \forall p\in {\mathrm{\Delta }}_{i,j},({\mathit{w}}^{\mathrm{a}}{)}_p=({\mathit{w}}^{\mathrm{f}}{)}_p+{\displaystyle \frac{({\mathit{x}}^{\mathrm{a}}{)}_{i,j}-({\mathit{x}}^{\mathrm{f}}{)}_{i,j}}{\text{card}\left({\mathrm{\Delta }}_{i,j}\right)}}.\end{array}$$

This approach was not favored in this first study, as it seemed more likely to generate negative weights more often.

Another alternative would be to generate new particles so that their weights sum up to the updated density, possibly fewer or more particles. This could be more technically challenging to implement and requires implementing the assimilation scheme directly inside the dispersion model loop. However, it could also have the advantage of conveniently increasing resolution where there are high concentrations of plastics.

7 Conclusions

This paper presents a simple yet readily effective method to assimilate observations of plastic concentration data into a Lagrangian dispersion model and its first implementation called the Ocean Plastic Assimilator (v0.2). We apply it in a controlled environment to assess its efficacy. We study the impact of observation errors on the prediction accuracy and changed some of the dispersion parameters (A and ϵ) to evaluate the impacts of model errors. Finally, we apply the method to a more realistic case with a real-world circulation field and find that the method still performs well. The encouraging results indicate that it is an excellent candidate to perform data assimilation with real-world data over ocean gyres.

Thus, the Ocean Plastic Assimilator will be further developed at The Ocean Cleanup to assimilate plastic concentration data from the oceans and improve our cleanup operations in oceanic gyres. This method will undergo more research to develop its features and assess its efficacy when using real-world observations. We expect it to be used to assess the cleanup operations of The Ocean Cleanup in real time.

The simplicity of the developed data assimilation method means that it should be easy to generalize to various other popular open-source Lagrangian frameworks such as OceanParcels (Delandmeter and van Sebille, 2019) or MITgcm (Campin et al., 2020). Porting the data assimilation procedure to the Julia language is also being envisioned whereby one could leverage the newly developed IndividualDisplacements.jl package to carry out Lagrangian simulations of plastic concentrations (Forget, 2021).