An analogues based forecasting system for Mediterranean marine litter concentration

In this work we explore the performance of a statistical forecasting system for marine litter concentration in the Mediterranean Sea. In particular, we assess the potential skills of a system based on the analogues method. The system uses a historical database of marine litter concentration simulated by a high resolution realistic model and is trained to identify 15 meteorological situations in the past that are similar to the forecasted ones. Then, the corresponding marine litter concentrations of the past analogue days are used to construct the marine litter concentration forecast. Due to the scarcity of observations, the forecasting system has been validated against a synthetic reality (i.e. the outputs from a marine litter modelling system). Different approaches have been tested to refine the system and the results show that 20 using integral definitions for the similarity function, based on the history of the meteorological situation, improves the system performance. We also find that the system accuracy depends on the domain of application being better for larger regions. Also, the method performs well to capture the spatial patterns but performs worse to capture the temporal variability, specially the extreme values. Despite the inherent limitations of using a synthetic reality to validate the system, 25 the results are promising and the approach has potential to become a suitable cost effective forecasting method for marine litter concentration.


Introduction
The ubiquity of the plastic waste pollution in seas and oceans worldwide raises great concern in the society and the scientific community, as it poses a significant environmental and 30 socioeconomic threat (UNEP, 2009). In consequence, the analysis of the impacts of marine litter pollution on the marine life and ecosystems has become a hot topic on marine research in recent years (Maximenko et al., 2019;Van Sebille et al., 2020;Lebreton and Andrady, 2019;Soto-Navarro et al., 2021). Marine litter particles accumulate both in shallow and deep waters, and particularly in enclosed basins such as the Mediterranean Sea (Soto-Navarro et 35 al., 2020;Cózar et al., 2015), where the observed concentrations are in the same range of those measured in the great plastic patches formed in the subtropical gyres of the open oceans (Cózar et al., 2015;Law et al., 2014;Van Sebille et al., 2015). Moreover, risk analyses have shown that marine organisms in the Mediterranean basin can be highly impacted by marine litter pollution (Compa et al., 2019;Soto-Navarro et al., 2021). The starting point to analyze those impacts and 40 to establish suitable mitigation strategies is to understand the spatial distribution and temporal evolution of the marine litter particles. Unfortunately, to carry on that analysis solely based on observations is not feasible. The large spatial and temporal heterogeneities of the field campaigns, along with the lack of standardized observational protocols, do not allow a synoptic representation of the marine litter distribution (see Maximenko et al. (2019) for a thorough analysis of the marine 45 litter observations problems and proposed improvements). For these reasons, numerical modeling emerges as a fundamental tool to achieve a synoptic description of marine litter dispersion patterns and as the base for the forecasting systems that would reproduce its spatial variability and time evolution.
Marine litter forecasting systems are usually based on the combination of two different numerical 50 models (Lebreton et al., 2012;Van Sebille et al., 2015;Maximenko et al., 2012). On the one hand, an ocean circulation forecasting system is implemented to provide ocean currents. On the other hand, a lagrangian model uses those currents to simulate the advection and diffusion of passive particles in the ocean that mimic the evolution of marine litter. In the Mediterranean, several studies using this methodology have been carried out using current fields from high resolution 55 regional models covering the whole basin (Liubartseva et al., 2018;Macias et al., 2019;Mansui et al., 2015;Soto-Navarro et al., 2020 ) or specific regions such as the Adriatic, the Tyrrhenian or the Aegean (Politikos et al., 2017;Fossi et al., 2017;Liubartseva et al., 2016;Palatinus et al., 2019). This modelling approach is considered to be the most accurate choice for marine litter forecasting (Van Sebille et al., 2020) provided the marine litter inputs are correctly prescribed 60 (Liubartseva et al. (2018) , Soto-Navarro et al. (2020)).
The downside of developing a forecasting system based on the direct modelling approach is that it involves a high technical complexity and computational cost. In order to overcome this limitations, it might be possible to develop a fast and light forecasting system based on statistical methods. One choice would be the so called Statistical Downscaling Methods (SDMs) which 65 relies on determining statistical relationships between large scale variables (usually atmospheric patterns) and local variables. They are broadly used in atmospheric modelling to forecast the evolution of local variables from large scale atmospheric models. The advantage of the SDMs is that the mathematical relationship derived by the model between the local and the large scale variables is valid not only for the present climate, but can also be used to estimate the future 70 evolution of the local variables. In summary, the SDMs provide a simplified 'static' methodology to forecast the evolution of local variables without the need of running a complex dynamical models. There are numerous downscaling methodologies based on different statistical properties.
Among them, the analogues method (Lorenz, 1969) is the most broadly used due to its simplicity and accuracy (Grouillet et al., 2016). This technique assumes that similar (or analogues) 75 atmospheric patterns over a given region, represented by large scale atmospheric variables or predictors, lead to similar local meteorological outcomes (or predictands) in a particular location.
This assumption provides a simple algorithm to downscale the local occurrence of the variable of interest from a given large scale atmospheric pattern (see section 2.1 for a detailed description).
In general, it has been shown that the analogues method performs as well as other more 80 sophisticated downscaling techniques (Zorita and von Storch, 1999), indicating that it is an efficient alternative for many downscaling problems. Its main advantages are that is nonparametric (i.e. no assumptions are made about the distribution of the variables used as predictors), non-linear (i.e. it can take into account the non-linearity of the relationships between predictors and predictands), and it is spatially coherent (i.e., preserves the spatial covariance 85 structure of the local variables). The analogues method has been satisfactorily applied in the Mediterranean region not only for the downscaling of meteorological or hydrological variables such as precipitation or river runoff (Grouillet et al., 2016;Wu et al., 2012;Caillouet et al., 2016), but also for the reconstruction of sea surface temperature in the glacial period (Hayes et al., 2005), the assimilation of satellite derived sea surface height (Lopez-Radcenco et al., 2019) and the 90 projection of complex climatic impact indices such as the fire weather index or the physiological equivalent temperature (Casanueva et al., 2014).
In this study, we explore the feasibility of a marine litter concentration forecasting system based on the analogues method. In particular, the surface marine litter concentration is linked to the atmospheric patterns during a reference period. Then, during the forecasting phase, the forecasted 95 atmospheric situation is compared to those realized during the reference period to identify analogue situations. The marine litter concentration during those analogues situations is considered to be a good approximation of the marine litter concentration that will occur during the forecasted date. As this is a new approach never tested before for marine litter dispersion, the first step has been to run several tests to fine-tune the methodology and to characterize its limits 100 of validity. Ideally, the tuning and validation of the method should had been done using in-situ observations but, unfortunately, the available marine litter concentration datasets are too scarce and this was not possible. Therefore, in this exploratory study, we have used numerically simulated marine litter concentration fields for the development and validation of the system. The rest of the paper is organized as follows. In section 2, the statistical method, the datasets used 105 and the different choices tested are introduced. In section 3, the model results are presented and discussed, and, finally, some conclusions about the capabilities of this new approach are outlined in section 4.

The analogues method 110
The implementation of the analogues method requires two sets of data. First, we need a reference dataset of the variables that describe the atmospheric patterns over the region of study, the so called predictors (X). The second reference dataset consists on spatial patterns of the variable of interest for the same period for which the predictors are available. In our case, those predictands (Y) would be the marine litter concentration fields. Once they are defined, the methodology is 115 based on the assumption that if two predictors are similar (X 1 ~ X 2 ) the corresponding predictands would also be similar (Y 1 ~Y 2 ). Then, in order to obtain a forecast of the marine litter concentration for a given date (Y fcst ), what we can do is to use instead a forecast of the predictor for the same date (X fcst ). In particular, we look for k analogue dates within the reference period (t an,k ) in which the predictor patterns are similar to the forecasted one (X (t an,k ) ≈ X fcst ). Then the value of the 120 variable of interest is estimated as a function of the predictands corresponding to the selected analogue dates Y fcst = f (Y(t an,k )). A scheme of the model algorithm is shown in figure 1a.
In our case, the predictors used to characterize the atmospheric conditions will be the Sea Level Pressure (SLP) and the wind speed (U 10 , V 10 ). These two variables have been successfully used to forecast ocean surface dynamics (Wang et al., 2010;Martínez-Asensio et al., 2016), so it is 125 reasonable to think that they may be also good to forecast marine litter concentration as it is mainly driven by ocean currents. The reference dataset for the atmospheric situation is obtained from an atmospheric reanalysis (see section 2.3). Regarding the reference dataset for the predictand, we use the marine litter concentration outputs from the modeling system developed by Soto-Navarro et al. (2020) and described in section 2.4. 130

Algorithm implementation
The first step to implement the analogues method is to define a cost function, JM, that measures the similarity between different meteorological situations. Then, for the forecast day (t fcst ) we estimate how close is the meteorological situation of that day to the rest of the days in the reference dataset by computing JM for the whole reference period. Those days with the lowest JM values 140 are selected as the analogue days ({t an,k }, see figure 1b for an example). For the definition of JM, the most popular choice is to use the Euclidean distance or root mean square error difference (RMSED) (Zorita et al., 1995;Cubasch et al., 1996;Gutiérrez et al., 2013), although other metrics based on different statistics can also be used. Here we have tested 4 different definitions for the cost function JM: 145 So, the similarity between meteorological situations is assessed either in terms of the sea level 150 pressure (SLP, JM 1 ), the 10-m winds (U 10 , V 10 ; JM 2 ), a normalized combination of both (JM 3 ) or the cumulated values of JM 3 during a period (Δt) before the reference day (JM 4 ). In our case, Δt has been set to 7 days. Note that the horizontal bars indicate spatial averages for (JM 1 ) and (JM 2 ), while <> in (JM 3 ) denotes temporal mean.
In a second step, we identify the analogue dates as those with the lowest values of JM. We keep 155 those dates in which JM is lower than the 1% percentile of all JM. Then, the marine litter concentration maps (C) obtained in the reference dataset for those days are combined to produce the forecast concentration map (C fcst ). In our case we use the median to reduce the influence of extreme concentration values close to marine litter sources:

Reanalysis data for the atmospheric fields
The period considered for the implementation of the analogues method is 2003 -2013, which coincides with the period simulated by the marine litter dispersion model (as described in the following section). The climatic dataset necessary for the model reference period is based on the ERA5 reanalysis dataset, available at the Copernicus Climate Change Service (C3S) web platform 165 (https://climate.copernicus.eu/climate-reanalysis). All the information regarding the ERA5 characteristics can be found on the C3S website.
Two variables have been considered for the characterization of atmospheric patterns forcing the marine litter dispersion: the wind speed at 10 meters height (U 10 , V 10 ) and the sea level pressure (SLP). Daily mean values of these variables over the Mediterranean Sea were downloaded and 170 processed for the whole period. The spatial resolution of the atmospheric data is 0.25º (~25 km) and cover the whole Mediterranean basin and the region of the North Atlantic adjacent to the Iberian Peninsula. Figure

Marine litter concentration data 180
The marine litter concentration data is obtained from the simulations performed by Soto-Navarro et al., (2020), as they are considered to be among the most realistic for the Mediterranean Sea.
Due to the relevance of the quality of the marine litter concentration data, some details on the modelling system are presented below and more information can be found in Soto-Navarro et al., (2020). 185 The system is based on two components, a regional high resolution circulation model (RCM) reproducing the 3D current velocity field in the Mediterranean (NEMOMED36), and a lagrangian model that simulates the evolution of floating particles (Ichthyop 3.3).
The hydrodynamical model used to simulate the Mediterranean current field is an implementation on the NEMO model, with a spatial resolution of 1/36 degrees (~ 3 km) in a domain that covers 190 the whole Mediterranean. The atmospheric forcing is a dynamical downscaling performed by the APEGE-Climate model using spectral nudging, namely ARPERA (Herrmann and Somot, 2008).
Note that the forcing of NEMOMED36 (ARPERA) is not the same that the one used to characterize the meteorological situations (ERA5). Although both datasets are very similar, they are not exactly the same, thus mimicking the inaccuracies that atmospheric forecasts will 195 inherently have.

Experiments 225
As mentioned before, there are no suitable observational datasets to validate the forecasting system. Homogenized datasets covering a long period of time would be required for this task.
Although there are some efforts to develop new databases (Maximenko et al., 2019), up to our knowledge, there are no such datasets in the Mediterranean yet. Thus, in order to have a first assessment of the quality of this methodology we have to use the concentration marine litter maps 230 from the database as a "virtual reality" and compare the forecast (C fcst ) with the C in the database for the forecast date (C(t fcst )). We are aware that this may produce overoptimistic results and this issue will be discussed below.
To define the forecast day, we pick any date from the reference period and forecast the marine litter for that day using all the data available except for a week before and after of the forecast day 235 to avoid spurious good results due to autocorrelation. This has been repeated for all the days in the reference period (3650) and several statistical metrics have been computed to assess the skills of the method.
To test if the model shows different skills depending on the domain of application, we have applied the method to seven different regions: the whole Mediterranean, the eastern and western 240 basins, and in the Gulf of Lions, the region around the Balearic Islands, the Adriatic Sea and the Aegean Sea (see figure 2). In each case, the analogue days have been defined using only data on the selected region.
Additionally, we have tested if the skill of the method depends on the time scales of the marine litter concentration variability. So, in addition to use the marine litter concentration dataset, we 245 have used two filtered versions of it, separating those processes above and below 15 days (Chi-freq and C lo-freq ).
Finally, for completeness, we propose three additional models for the forecasting. First, we forecast the concentration change in 7 days (Δ 7d C). The underlying idea is that the meteorological situation could be a better predictor of the rate of change than of the absolute value (e.g. winds 250 may determine the changes in the concentration rather than the absolute value). The second one is to simply assume 7-days persistence as the forecasting model (I.e In summary, we have tested 4 configurations of the model over 7 different regions to forecast C , C hi-freq and C lo-freq 260

Quality assessment
Several diagnostics are used to characterize the quality of the forecasts in the different experiments. The first one is the root median square error (RMEDSE): We have chosen this parameter instead of the root mean square error to reduce the overall impact 265 of outliers linked to very high concentration values close to marine litter sources. Complementary we also compute the temporal correlation ρ: where Cov represents the covariance and the standard deviation. Additionally, we compute the RMEDSE ratio (RR) which is defined as the ratio between the RMEDSE of the forecast (eq. 6) 270 and the RMEDSE computed using all the days in the database, RMEDALL: The lower the value of RR is, the better the forecast is. Values of RR close to 1 means that the quality would be the same than using any random day, so the forecast is not providing any new information. RMEDSE, ρ and RR are computed spatially and/or temporally. 275

Time variability
The temporal correlation and the RR of the marine litter concentration reconstruction using Concerning the different cost functions used to identify the analogue situations, the performances using only SLP (JM 1 ) or only wind (JM 2 ) are very similar. Using both variables the quality slightly increases (JM 3 ) and becomes significantly better when using the 7-days average (JM 4 ).
For model 1 (forecasting the concentration), the averaged correlation using each cost function is 0.24, 0.25, 0.28 and 0.35 while the averaged RR is 0.93, 0.93, 0.90 and 0.86, respectively. The 290 forecasting of the concentration change is worse for all cost functions, with averaged correlation values ranging from 0.08 to 0.19 and RR ranging from 1.00 to 0.98. In the light of these results, from now on, we will only consider the results of the analogues-based forecast models that use the cost function JM 4 (i.e. the one considering the 7-day averaged differences). Using it for forecasting the marine litter concentration we obtain correlation values ranging from 0.20 to up 295 to 0.60 depending on the region. When forecasting the marine litter concentration change the values range from non-significant to 0.40 (see Figure 4).  For completeness, we also include an example of the concentration time series for the reference 315 and models 1, 3 and 4 for a point where the forecasts perform well (Figure 6a). It can be seen that Model 1 is well correlated with the reference, showing a good chronology of events although being unable of capturing the concentration peaks. During those periods, the analogues-based forecast largely underestimates the reference values. Models 3 and 4 show almost identical good results, as far as persistence is enough to capture most of the variability. The underlying reason 320 for this success is that, at this location, the changes of marine litter concentration are relatively slower, so assuming persistence can be a good predictor. For comparison, the time series for a point where the models perform poorly are shown in Figure 6b. In this case, the analogues-based forecast is unable to capture any variability and it basically produces the mean value. The other two models are able to follow the variability, although in this case the skills are lower than in the 325 previous case. The reason is that in this point the marine litter concentration varies more rapidly, so assuming the persistence is not as good predictor as it was in the previous location.

Spatial Variability
A complementary view of the performance of the different forecasting models can be obtained looking at the marine litter concentration anomalies (i.e. with respect to the temporal mean) in 335 given dates. In Figure 7, we show the results for a date when models show good agreement with the reference (spatial correlation values are 0.70, 0.76 and 0.78 for models 1, 3 and 4, respectively). All three models are able to identify the areas of high and low concentration.
Maximum values in the north of the Balearic Islands, the Gulf of Gabes and south of Italy and minimum values in the Adriatic Sea, the Algerian basin and the easternmost part of the Med are 340 well captured. The analogues-based forecast (Model 1) shows smoother patterns with less low extremes. This is in good agreement with what has been seen in the time series in Figure 6, suggesting that this model has difficulties to capture very high concentration values. Regarding the persistence-based models, for this particular date, they perform very well capturing not only the large scale patterns but also the local features. Looking at a date when the performance is 345 lower something interesting appears. Although the spatial correlation of Model 1 is not significant (Figure 8b), the large scale features seem to be well captured. However, the small scale features are clearly not captured which degrades the spatial correlation. This would also support the previous finding reinforcing the idea that the analogues-based forecast performs better for the large scale features. In places or dates where/when the small scale features become dominant, the 350 performance of the model drops. The time series of the spatial correlation and spatial RR at each time step are presented in Figure   9. The results are similar for the three models forecasting the marine litter concentration (Models 1, 3 and 4). The skills of the forecasts show a high temporal variability with correlation values 360 ranging from 0.5 to almost 1 and an averaged value of 0.78, 0.81 and 0.84, respectively. For RR the values range from 0.3 to more than 1 with an average value of 0.76, 0.79 and 0.71, respectively. This diagnostic also confirms that the best model is the one combining the persistence with the forecast of the concentration change.

Regional dependence of the forecasting skills
The methodology has also been applied to different domains. That is, the cost function, JM, has 370 been computed in the regions defined in figure 2 and the validation has been performed looking only at the marine litter concentration in those regions. In general, better results are obtained when the analogues-based forecasts are applied to a larger region (see Table 1 and Table 2  The spatial diagnostics have also been computed applying the models to different domains (Table   3). In this case, the analogues-based forecast of concentration (M1) show average spatial 400 correlations higher than 0.62 when applied to any region reaching up to 0.94 in the Aegean Sea.  Table 3. Temporally averaged regional correlation and RR of the different forecasting models (M1-M4) applied in different regions (see Figure 2).

410
It is worth mentioning that we have also tested other options for the cost function like using different temporal averages or using correlation as similarity metrics but no significant differences have been found. Also, we have tried to change the criterion to define the analogue days. Instead of identifying as analogues those days with JM lower than the 1% percentile of the whole JM time series, we have used less restrictive criteria (5% or 10%). In both cases the results worsened. 415

Discussion and Conclusions
The analogues-based forecasting technique has been applied to marine litter concentration for the first time, up to our knowledge. It has proven to be very inexpensive and relatively easy to setup, so it is an alternative to direct modelling worth to be considered. A key step in the set-up is to select a suitable cost function and the best threshold to identify the analogue meteorological 420 situations. In our case, it seems that using integral definitions for the cost function improve the results. In other words, it is better to identify the analogue days based on the history of the meteorological situation. Probably, using a different averaging time for each domain would allow increasing the skills of the analogues-based model. However, this fine tuning is out of the scope of this paper, as far as there are no suitable observations to validate it, as it will be discussed later. 425 The quality of the analogues-based forecasts depends on the region of application. Our results suggest that the larger the region of application the better, as we get better results for the whole Mediterranean or for the East/West basins than in smaller local areas. A hypothesis for explaining this result is that using the atmospheric situation as a predictor may not be suitable to capture small scale features (e.g. those related to ocean currents or the interaction with coastlines). Further 430 tests including other predictors could be done to refine the method, including ocean currents, for instance.
Another important point is that the method struggles to capture the extreme values as it produces smooth spatio-temporal patterns of marine litter concentration. Therefore, in locations or regions where short intense events or small scale features dominate the variability, the method performs 435 worse. This is also one of the reasons why the temporal skills are relatively low (i.e. temporal correlation and RR, see section 3.1). Conversely, if instead of the time variability, what are aimed at are the spatial structures, the method shows high skills being able to locate relative maximum and minimum (see section 3.2).
We have also shown that persistence is a very good predictor almost everywhere. This is because 440 the marine litter concentration changes relatively slowly (i.e. the system has a several days memory), at least at the spatial scales solved by the reference dataset. This means that if reliable information was available (e.g. from a monitoring program), this could be used as a first guess of the marine litter concentration several days later. Complementary, the analogues-based method has also been applied to forecast concentration change. In this case the results were significantly 445 poorer both to capture the time and spatial variability. However, it can be useful to improve the persistence-based forecasts.
Regarding the reliability of the analogues based forecasts that could be generated from this reference dataset, its quality would directly depend on the accuracy of the reference dataset. In our case this dataset comes from the outputs of a realistic modelling (Soto-Navarro et al., 2020). 450 However, the model may have some shortcomings as its spatial resolution, beaching parameterization or realism of marine litter sources. Consequently, the forecasts would be, in the best case, as good as the model outputs are. Therefore, it would have been better to validate the different forecasting models against actual observations. Unfortunately, the lack of observations with a suitable spatial and temporal coverage prevents from doing it. In the future, it would be 455 worth setting up a monitoring program with enough spatial and temporal resolution that would allow generating a comprehensive enough reference dataset. This dataset could be used to train the analogues-based forecasting system and to validate other existing systems.
In any case, it is worth noting that the validation of the methodology can be considered as robust.
For that purpose, it is not required that the reference dataset is an accurate representation of the 460 actual marine litter concentration. Only the statistics of the marine litter concentration spatiotemporal evolution has to be reproduced. And in that sense, the model integrates the effects of a realistic atmospheric forcing and a realistic ocean current field. So, it is expected that the statistics of the marine litter concentration field is realistic enough. This extent should also be confirmed by a comprehensive observational dataset, at least in certain regions. 465 In conclusion, the analogues-based model presented here has potential to become a suitable cost effective forecasting method for marine litter concentration. It could be easily implemented in any region of the world where a realistic reference dataset is available. In those regions where the large scale marine litter concentration patterns dominate the variability the method will probably work better than in regions where the variability is dominated by small scale structures. 470

Code and data availability
The code and data required to implement the model described in the paper and to reproduce the results can be publicly accessed at Jordà and Soto-Navarro, (2022). Additionally, the atmospheric fields can be downloaded from the Copernicus portal (https://climate.copernicus.eu/climatereanalysis). 475

Author contribution
Both authors (GJ and JSN) have contributed equally to the design of the study, the coding of the modelling system, the performance of the simulations, the analysis of the results and the preparation and revision of the manuscript.

Competing interest 480
The authors declare that they have no conflict of interest.