Strategies for simulating the drift of marine debris

Modelling the drift of marine debris in quasi-real time can be of societal relevance. One pertinent example is Malaysia Airlines ﬂ ight MH370. The aircraft is assumed to have crashed in the Indian Ocean, leaving ﬂ oating wreckage to drift on the surface. Some of these items were recovered around the western Indian Ocean. We use ocean currents simulated by an operational ocean model in conjunction with surface Stokes drift to determine the possible paths taken by the debris. We consider: (1) How important is the in ﬂ uence of surface waves on the drift? (2) What are the relative bene ﬁ ts of forward- and backward-tracking in time? (3) Does including information from more items re ﬁ ne the most probable crash-site region? Our results highlight a critical contribution of Stokes drift and emphasise the need to know precisely the buoyancy characteristics of the items. The di ﬀ erences between the tracking approaches provide a measure of uncertainty which can be minimised by simulating a su ﬃ ciently large number of virtual debris. Given the uncertainties associated with the timings of the debris sightings, we show that at least 5 items are required to achieve an optimal most probable crash-site region. The results have implications for other drift simulation applications.


Introduction
Drift observations of objects near the surface of the ocean have led to fundamental advances in the understanding of ocean dynamics. Most notably, Nansen's observations of ice drift in the Arctic led to the formulation of Ekman theory, which plays a central role in understanding the dynamics of ocean gyres, fronts, and upwelling, among others. Since the latter part of the twentieth century, satellite-tracked drifting buoys continuously provide a wealth of information on ocean circulation and in-situ properties (Lumpkin et al. 2017).
Supplementing these observations, simulations of pathways of floating objects in the marine environment are routinely undertaken, in particular for assessing the dispersion and accumulation of marine litter and pollution (e.g. Hardesty et al. 2017;Robinson et al. 2017). Such assessments require model data that allow for the connection between the coastal environment (as primary source of anthropogenic litter or pollution) and the open-ocean, and thus span wide spatial and temporal scales. The model data generally stem from simulations that, to a large extent, faithfully represent the large-scale upper ocean circulation of the past decades. However, given the high spatial and temporal variability of real-world surface flows resulting from features such as mesoscale eddies, the results of such studies must be interpreted in a probabilistic context, often with large uncertainties.
Here, we provide insights on surface drift simulations of objects for quasi-real-time applications. Instead of employing simulated data from the past decades using ocean models driven by reanalysed surface fluxes, we utilise an up-to-date description of simulated surface velocities which includes the assimilation of direct observational data. As a case study, we consider debris of the Malaysian Airlines (flight MH370) aircraft that were collected in the western Indian Ocean over the period July 2015 -June 2016.
On 8 March 2014, flight MH370 disappeared on its way from Kuala Lumpur, Malaysia to Beijing, China with 239 people on board. Initial analyses from radar and other communications pointed to a deviation of the planned flight path towards a southern route in the Indian Ocean (ATSB 2017). Detailed analysis of satellite communications, provided in the form of handshakes between the aircraft's engines and satellites, indicated that the plane had lost contact along the '7th arc' around the position of the satellite, which extended from Java, Indonesia, to the southern Indian Ocean, southwest of Australia (Figure 1, ATSB 2017). However, the precise location of the aircraft's last position was not known. In the following months, the Joint Agency Coordination Centre, led by the Australian government, started the search for the aircraft around the 7th arc. The unsuccessful search was halted in January 2017 and spanned an area in excess of ∼120000 km 2 (Figure 1).
Several items of debris were recovered around the Indian Ocean and attributed to the MH370 aircraft ( Figure 1). The first and most significant piece was a flaperon found on La Réunion on 29 July 2015, about 16 months after the disappearance of MH370. Presuming that this item had been drifting with the ocean currents after the crash, could its discovery help the search for the missing plane? Or more generally, can 'backtracking' an object's path at the ocean surface (i.e. following it backwards in time from where it was found) provide meaningful insights on its origin?
In answering this question, several issues arise. We outline here three key considerations that we examined through the MH370 case study.
(1) The drift of any floating object on the ocean surface is influenced by: (1) surface ocean currents, (2) Stokes drift (nonlinear net drift near the ocean surface which results from surface waves), and (3) windage resulting from direct wind drag (tangential shear stress), form drag (pressure differences across the item), or the 'sail' effect (resulting a force at an oblique angle to the wind direction). The geometry of the object determines the respective contribution of these factors to the total drift. (2) Describing motion from a particle-following Lagrangian perspective allows for both forwardand back-tracking in time. However, these two approaches do not necessarily yield the same result. We compare both approaches, even though for the flaperon, a back-tracking approach may seem more intuitive since presuming any crash site a priori might not be sensible. Figure 1. Geographical distribution of the discovery locations of 9 items of debris ordered according to their discovery dates. The colour scale shows the ocean depth (m) for the region over which virtual objects were released (see text for details). The maximum range is a perimeter delineating the extent the aircraft could have flown given its fuel load and consumption. The 7th arc is the possible position of the aircraft during its last successful automatic satellite communication. The hatched patch roughly shows the area of the underwater search undertaken by the Australian Transport Safety Bureau (ATSB) between April 2015 and January 2017.
(3) Hypothetically, a more robust probable crash site could be obtained by including more items of debris into the analysis. We assess this hypothesis under the consideration that for most items, their date of discovery might differ substantially from the true date of arrival at their respective locations.
While we focus on the MH370 case study to illustrate our analyses, our results are valid for other applications of object drift simulations in the ocean.

Materials and methods
Over the period July 2015 -June 2016, 20 pieces of debris were collected along the western Indian Ocean coasts, of which 9 were identified as 'confirmed' or 'almost certain' to belong to the MH370 aircraft ( Figure 1, MOT 2018). So as to properly simulate their possible drifts, we combined surface currents obtained from the Copernicus Marine Environment Monitoring Service (CMEMS) operational model with matching simulated Stokes drift. We first describe the model products we employed and then elaborate on the Lagrangian methodology.

Model data
We used daily-mean surface velocities for the period 1 March 2014 to 26 July 2016 from the CMEMS global ocean 1/12°physics analysis and forecast model in its most accurate delayed mode. The system uses a global ocean model based on the NEMO code (v3.1; Madec 2012) in the ORCA12 configuration. It covers the global ocean at 1/12°nominal resolution (∼9 km grid size in the Indian Ocean) and has 50 levels in the vertical (surface layer 1 m thick). The model is forced by 3-hourly European Centre for Medium-Range Weather Forecasts (ECMWF) operational winds and corresponding heat and freshwater fluxes. It captures all effects of these fluxes in driving the ocean currents, including the Ekman drift from the wind, but excludes wave effects. The model is run operationally, assimilating observational data. Using a Kalman filter with the SEEK formulation (SAM2v1) and bias correction (3D-Var) with incremental analysis updates, a range of satellite and in-situ data are used to update and correct the simulated ocean state (Lellouche et al. 2013(Lellouche et al. , 2018. In particular, the use of sea surface temperature and sea surface height from satellites combined with the mean dynamic topography provide a more accurate description of the upper-ocean hydrography and velocity, resulting in a state-of-the-art description of the ocean surface properties at any given time. In the Indian Ocean the sea level anomaly bias is less than 0.005 m 1 and the mean error in surface (15 m) velocity, assessed against independent observations, is ∼0.1 m/s (Lellouche et al. 2018). The synoptic evolution of mesoscale features such as eddies is also more realistically reproduced by the operational model than by a non-assimilative hindcast.
Stokes drift was taken from the ECMWF High RESolution WAve Model (HRES-WAM; Breivik et al. 2016;ECMWF 2017). Using an atmospheric forcing consistent with that used for the ocean model described above, the wave model explicitly simulates Stokes drift at 1/4°resolution. We interpolated daily-mean Stokes drift velocities onto the ORCA12 grid and added these to the ocean model surface current velocities.

Lagrangian tracking of objects
Pathways of objects passively transported by a given flow can be described from a Lagrangian perspective. In this frame of reference, an object starting at a given position and time is moved from that position either forward or backward in time using the Eulerian velocity description of the flow field for each time step. Lagrangian simulations of pathways within ocean models are now common practice (van Sebille et al. 2018) and have been previously performed to trace pathways of water masses (e.g. Durgadoo et al. 2013Durgadoo et al. , 2017Rühs et al. 2013) or marine litter/pollution (e.g. Robinson et al. 2017).
The Lagrangian experiments for this study were performed using the Ariane software (v2.2.6, Blanke and Raynaud 1997). Ariane simulates the advection of objects within a simulated volume-conserving velocity field by displacing them along analytically computed streamlines. It carefully considers the C-grid layout of NEMO, and has already been successfully run using velocities simulated by the ORCA12 configuration (e.g. Durgadoo et al. 2017;Robinson et al. 2017).
For the MH370 case study, the following assumptions were made: . All items of debris were simulated as dimensionless virtual objects passively floating on the ocean surface (at a constant depth of 0.5 m), in other words being incapable of self-propulsion, and experiencing zero drag. . Windage on all items of debris was assumed negligible.
For the flaperon, buoyancy analyses and the presence of barnacles suggests that it floated in a near-horizontal, slightly submerged position and thus with negligible windage (Daniel 2016; ATSB 2017). . By using daily-mean velocity fields the stochastic effects of quickly changing weather events on the drift of objects is neglected. . It is not known how long each item of debris had been on the beach before its discovery, or how long it had been drifting in coastal waters before being washed ashore. Additionally, CMEMs may not simulate the full spectrum of coastal processes faithfully. To account for such spatial and temporal uncertainties, we considered virtual objects within a 1°x 1°box around each location and covering a 30-day 'leeway' period prior to their discovery (Table 1). . The authorities assumed that the aircraft crashed somewhere in the south-eastern Indian Ocean. Here, we first considered the entire Indian Ocean basin as the possible crash site. Subsequent refinements were undertaken based on the maximum flight range of the aircraft or along a swath around the 7th arc.
Two reference experiments were undertaken: 1. Back-tracking: Virtual objects were initiated uniformly (1/12°spacing) within a 1°× 1°box around La Réunion at hourly intervals for 30 days prior to the discovery date (∼112000 independent objects). Sensitivity tests were performed to examine the effect of including various proportions of the full Stokes drift, different leeway periods, and the number of particles released. These are detailed in Table 1.

Subsampling and probability calculations
For each back-tracking experiment, starting at La Réunion, the positions on 8 March 2014 of all virtual objects within the domain were recorded. These were binned and counted on a 1.5°× 1.5°regular grid. A probability distribution was achieved by dividing the count in each bin by the total number of virtual objects.
For each of the forward-tracking experiments the following procedure was undertaken: (1) For each debris discovery location X i , 1 ≤ i ≤ 9 (a) Isolate those virtual objects that, within the prescribed leeway period, reach a 1°x 1°box around X i . (b) The start positions (on 8 March 2014) of the subset of objects in (a) are binned and counted on a 1.5°x 1.5°regular grid. (c) A probability distribution P(X i ) is calculated by dividing the count in each bin by the total number of virtual objects isolated in (a).
(2) A joint probability distribution is calculated as P(X 1 AND X 2 … AND X 9 ) = P(X 1 ) x P(X 2 ) x … x P(X 9 ).
The joint probability yields a distribution within the overlapping region of all individual probability distributions. We normalised the joint probability distribution by dividing by ∑[P(X 1 ) x P(X 2 ) x … x P(X 9 )].
We portray the (joint) probability distributions by delineating a 'most probable' region which contains the highest probabilities and which give an accumulated probability of 0.75 (i.e. the region we assess as most likely to contain the crash site).
The above-mentioned strategy highlights two ways to account for uncertainties that result from the chaotic nature of the (simulated) surface ocean flow. Firstly, by Table 1. Details of the experiments performed. The leeway (given as days prior to the reported discovery date) accounts for not knowing how long the items of debris have been ashore or drifting in coastal waters before their discovery.
Release: uniform within the Indian Ocean basin using a relatively coarse 1.5°× 1.5°binning, corresponding to an area roughly ranging between 18000 and 28000 km 2 (depending on the location of the bin). Secondly, by simulating the trajectories of a sufficiently large number of objects (assessed through the sensitivity tests which were performed, Table 1). The resulting probability distributions provide a quantification of the associated uncertainty.

Results and discussion
Contribution of wave-induced motions to surface drift The surface flow within the Indian Ocean can be broadly separated into three regions. The equatorial and northern Indian Ocean is characterised by seasonal reversals of winds and currents (Schott et al. 2009). The southern Indian Ocean consists of a broad sluggish flow associated with the subtropical gyre (New et al. 2007) and two fast flowing western boundary currents off the east coasts of Madagascar and South Africa (Lutjeharms 2006). Finally, south of the Agulhas Return Current, which lies near 42°S, the Indian Ocean sector of the Southern Ocean is strongly influenced by persistent westerly winds (Durgadoo et al. 2008). This diversity in regimes results in a wide range of instantaneous (daily-mean) simulated surface currents (Figure 2a). The simulated surface Stokes drift for 8 March 2014 resembles the pattern of the prevailing winds ( Figure  2b). A non-negligible contribution O(0.2 m/s) from these wave-induced velocities is seen in the southern and eastern Indian Ocean, adding an anticlockwise component to the combined velocities between approximately 10-40°S which is northwestward in the region west of Australia (Figure 2c). On average and for regions outside the subtropical gyre, Stokes drift accounts for ∼40% of the combined velocities (Figure 2d). Within the subtropical gyre, Stokes drift and surface ocean currents are comparable.
We assess the impact of Stokes drift on the drift of objects by using the back-tracking experiments from La Réunion ( Table 1). The probability distributions of the positions of the objects on 8 March 2014 under different percentages of Stokes drift are shown in Figure  3. The overall anticlockwise pattern of the wave-induced velocities is reflected in the progressive shift of the distributions from the northern to the southwestern Indian Ocean, as the percentage of the simulated surface Stokes drift increases from 0 to 150%. The most probable area increases ∼4 fold in size to ∼2 million km 2 when accounting for Stokes drift at 100%. The ± 50% sensitivity experiments for Stokes drift provide a measure of the uncertainty in exactly how waves interact with the debris (in this case the flaperon). For example, it is uncertain how shorter waves impact the debris and over what effective depth Stokes drift should be applied. Indeed, a substantial percentage of the surface Stokes drift, typically about a third (Breivik et al. 2016), results from the high frequency tail of the surface wave spectrum. This component decays rapidly with depth, with a decay scale that may be as short as 0.1-0.2 m. It is thus questionable whether an item of debris 'feels' the full surface Stokes drift or only part of it. Other factors such as the reflection of short high-frequency waves from the debris could enhance the drift effect (Longuet-Higgins 1977). Additionally, some laboratory work (e.g. Huang et al. 2011) has suggested that waves may cause objects to move faster than Stokes drift.
Windage is ignored in this study. Daniel (2016) analysed the floating characteristics of the flaperon and concluded that (without barnacles) such 'windage' could be ∼3% but that with the presence of barnacles, as actually observed, the flaperon would float slightly submerged and with 0% windage. Given that barnacles take a finite time to grow, this is an uncertainty in our study. Additional windage would in general tend to push the most likely region further back around the gyre, given that the prevailing winds are eastward near 40°S and then northwestward in the eastern regions of the gyre. Given that Stokes drift is typically 1-2% of the wind speed (Smith 2006;Ardhuin et al. 2009), an additional 50% in the Stokes term (e.g. from 100% to 150%) would be comparable to an additional windage of 0.5-1.0% of the wind speed. However, the (direct) effect of wind on surface drift might not be aligned to that of waves. Therefore, the comparison between Figure 3c and d only gives a plausible difference in distribution should such an additional windage be included.
Noting the above-mentioned uncertainties, we choose to proceed with the advection of virtual debris with the full surface Stokes drift (100%) in the following as being the most realistic case.

Forward-vs back-tracking
The above back-tracking Lagrangian experiment is intuitive since it makes no a priori assumptions on the possible origins of the debris. Based on the knowledge of the Figure 3. Probability distribution on a 1.5°× 1.5°regular grid of object locations on 8 March 2014 from the La Réunion back-tracking experiment. For these, simulated surface currents combined with 0, 50, 100 and 150% of simulated Stokes drift (a-d respectively) were used. In colour, a most probable area is defined whose accumulated probability is 0.75, which is enclosed within the grey area (accumulated probability = 1). possible position of the aircraft during its last successful satellite communication (around the 7th arc), a forwardtracking experiment could also be devised with virtual flaperons starting around the arc. These two experiments would, however, not necessarily yield the same results.
Starting from the forward-tracking reference experiment (Table 1) where virtual objects were spatially uniformly released within the entire Indian Ocean on 8 March 2014, we isolated those objects that reached a 1°× 1°box around La Réunion between 30 June 2015 and 29 July 2015 (30 days leeway). This subset consisted of ∼27000 objects. Figure 4a shows the probability distribution of the start positions of these objects. The positions and timings of these subsampled objects around La Réunion were noted and subsequently used to initiate a reverse/back-tracking experiment to determine their respective positions on 8 March 2014. This experiment is similar to the back-tracking reference experiment listed in Table 1, but with fewer objects and non-uniform starting positions around La Réunion. The resulting probability distribution is shown in Figure 4b. That Figures 4b and 3c yield similar results (spatial correlation of 0.9 at the 95% confidence interval, c.i.) suggests that the number of virtual objects used for the reference back-tracking experiment is sufficient.
Quantitatively, the distributions shown in Figure 4 are not similar (spatial correlation of 0.5, c.i. = 95%; a Kolmogorov-Smirnov test rejects the null-hypothesis that the two distributions stem from the same underlying continuous population, c.i. = 95%, p∼10 −4 ). The Figure 4. Probability distribution on a 1.5°× 1.5°regular grid of object locations on 8 March 2014 from the forward-tracking to La Réunion (a) and the La Réunion reverse/back-tracking experiment (b). For these, simulated surface currents combined with 100% of simulated Stokes drift were used. In colour, a most probable area is defined whose accumulated probability is 0.75, which is enclosed within the grey area (accumulated probability of 1).
forward-tracking experiment (Figure 4a) produces a most probable area ∼2.5 million km 2 larger than that produced by the reverse/back-tracking experiments (Figures 4b and 3c). In theory, the analytical solution of the trajectory equation, as employed by the Ariane software, is unique for 3-dimensional non-divergent flows. Thus, integrating a trajectory forward in time and then using its final position and integrating a trajectory backward for the same period of time should yield identical results. In practice, however, tests show that full reversibility is compromised by slight numerical errors. Analytical advection in 2-dimensions within Ariane is achieved by setting the vertical velocity component of the flow to zero, hence forcing each object at each time step to remain at the given constant depth or density. Thus, a further small error is introduced at each time step which also accumulates as trajectories are calculated. The overall effect is that the solution is no longer unique and forward and backward trajectories differ. These factors are responsible for the differences between Figure 4(a) and (b). Other explicit time stepping methods for advection in 2-dimensions would yield similar results (van Sebille et al. 2018). Additional sources of error include slight inconsistencies in timing when initiating the reverse tracking. Such inconsistencies can be minimised by using velocity data stored at high temporal resolution and subsequently reducing the trajectory integration time step. Nonetheless, the qualitative similarity between the probability distributions resulting from both strategies shows that the net effect of the above errors is relatively small. This is achieved by having a sufficiently large number of independent virtual objects.

Inclusion of multiple debris discovery locations
To explore whether including additional items of debris can further reduce the uncertainties, we could either run 8 further separate back-tracking experiments or use the single forward-tracking experiment for all debris locations. The latter strategy ensures that each debris item is given equal weight, and is more consistent with the prognostic input fields. This is therefore the strategy we use here.
Following the two-step procedure described above, for each of the 9 debris discovery locations, individual probability distributions from the reference forwardtracking experiment were calculated (similar to that shown in Figure 4a). Their normalised joint probability distribution is shown in Figure 5a. The most probable region (with accumulated probability of 0.75) roughly spans an area of ∼1.0 million km 2 in the southwest Indian Ocean between 30-80°E and 35-40°S. Compared to using only the flaperon (item #1), this is a reduction in area of ∼90% (Figure 5b). The area of the overall distribution (with accumulated probability of 1) is reduced by ∼50% when debris from all 9 locations are considered. Figure 5b further shows a quasi-asymptotic behaviour after the 5 th item of debris has been considered; that is the area reduction beyond that point is relatively small. Here, the debris locations have been considered in the temporal order in which the respective items were reported. However, it is worth noting that this quasiasymptotic behaviour does not depend on the order in which the locations are considered (not shown).
As indicated previously, a 30-day leeway margin was chosen to account for the temporal uncertainty associated with not knowing precisely when each item of debris reached its respective location. We tested this assumption by altering the margin to 10 and 50 days. The resulting areas of the most probable region were found to be similar to that of the reference experiment (not shown), and to show a similar exponential decrease in area with increasing number of debris locations considered (Figure 5b, blue and red curves respectively, compared with the black curve for the reference experiment). These two sensitivity experiments provide an error bar to the reference experiment.
In this context, we consider what might cause the differences between the most likely region in Figure 4a (using only the flaperon) and that in Figure 5a (using all 9 items). Firstly, as each subsequent item of debris is included in the analysis, the most probable regions show a regular progression between Figures 4a and 5a (not shown), retreating from the central and eastern portion of the South Indian Ocean between 30-40°S and moving steadily to occupy a more westward position. Also, there is a general anticlockwise circulation in this region ( Figure  2d; New et al. 2007) which would bring surface objects from the west to the east (between 30-40°S) and thereafter northwards west of Australia, subsequently reaching the westward South Equatorial Current between 10-15°S which would be expected (on average) to bring the objects towards Reunion, then Madagascar and then to the African coast. Since the inclusion of further items of debris tends to drive the most probable region further backwards around this circulation pattern, we speculate that there could be significant delays between beaching and discovery for some of the later items, that would have the spurious effect of giving more time for the items to circulate and result in starting positions further back around the gyre than they should be. This speculation is supported by the observation that some items of debris appear to be discovered significanly 'out of sequence' with e.g. items 5 and 7 being found 8 and 10 months after item 1, even though they are 'upstream' and would usually be expected to be found earlier.
Finally, are ∼22 million starting objects used for the reference forward-tracking experiment sufficient? Reducing the number of starting objects by half yields probability distributions (spatial correlations of ∼0.8, c.i. = 95%, not shown) that fall well within the error bar that we defined above (Figure 5b, green curve), suggesting that the number of objects used for the reference forward-tracking experiment is sufficient.

Implications for the search for MH370
Coming back to whether the discovery of debris could assist the search for the missing MH370 aircraft, 2 additional clues might help to restrict the geographical area: (1) the maximum range the aircraft could have flown given its fuel load and consumption rate and (2) the possible position of the aircraft during its last Figure 5. (a) Normalised joint probability distribution on a 1.5°× 1.5°regular grid derived from all 9 debris locations on 8 March 2014 from the forward-tracking experiment, which uses simulated surface currents combined with 100% of simulated Stokes drift. A most probable area is defined whose accumulated probability is 0.75 (blue area), which is enclosed within the grey area (accumulated probability = 1). (b) The area (km 2 , left axis) and % area decrease from only considering the first debris location (%, right axis, lines with brown circles), as a function of the number of debris locations considered (taken in order as shown in Figure 1). All curves describe the region with an accumulated probability of 0.75, except for the thick dashed curve which shows the area of the entire distribution (accumulated probability of 1).
successful satellite communication (along the 7th arc). From Figure 5a it is clear that the region within the maximum range perimeter (and indeed around the arc) is associated with low probabilities, that is, the chance that all 9 items of debris originated from there is low. However, Figure 6 shows the result of refining the calculations for Figure 5a by only considering those trajectories which originated either within the maximum range (Figure 6a) or both within the maximum range and near the arc (Figure 6b, for which we used a ∼550 km swath on either side of the arc). In the latter case, the most probable region for the crash site around the Arc lies between 30-35°S and 80-95°E. We note in particular that the area searched by the Australian Transport Safety Bureau (here also binned onto a 1.5°× 1.5°grid) overlaps with the region of highest probabilities that intersects the arc (Figure 6b). At the time of writing, metadata was not available from the 2018 Ocean Infinity search, which extends to 25°S.
In the full list of debris recovered around the western Indian Ocean, 2 items were identified as 'likely' to originate from MH370, 8 as 'very likely', and 7 as 'not identifiable' (MOT 2018). None of these were considered in the present study. From Figure 5b, we deduce that the inclusion of more debris discovery locations in the Figure 6. Refinement of Figure 5a by considering (a) the maximum possible range of the MH370 aircraft and (b) a 5°swath along either side of the 7th arc between the island of Java and 45°S. For each panel, in colour, a most probable area is defined whose accumulated probability is 0.75, which is enclosed within the grey area (accumulated probability = 1). analysis is not likely to lead to a significant refinement of the most probable area. Also, with a 30-day leeway, item #8 (Antsikara) is representative of 7 other items retrieved around the same site in June 2016.
When considering the implication of our results for the search for MH370, we must acknowledge several sources of uncertainties, which led to the assumptions made in this study. The geometry of a debris determines its buoyancy characteristics (e.g. drift depth and windage) and hence influence its drift. These characteristics would be different for different pieces of debris. Here we demonstrated the usefulness of the methodology that we employed in including multiple items of debris. Our results could be refined further if, for all debris, their buoyancy characteristics and their time of arrival on land were known precisely.

Summary and conclusions
We undertook a detailed study of trajectory calculations in the Indian Ocean in connection with the disappearance of flight MH370 in 2014. We used a state-of-the-art ocean circulation model which includes the assimilation of available data to provide the best possible estimates of the ocean surface currents during the relevant period.
We showed the importance of including Stokes drift into the calculations. This is a nonlinear effect which results from the wave field in the ocean, and critically affects the results since this term can be of the same order as the wind-and buoyancy-driven surface ocean currents in large regions of the subtropical gyre. Entirely ignoring Stokes drift in the calculations can result in large errors, which in our case study clearly demonstrates. Although there are still uncertainties arising from the windage on the floating items of debris, we believe that this effect is relatively small (at least for debris item #1) and likely to be much less than the differences between using 50% or 150% of the Stokes drift term. It is our opinion that the scenarios including Stokes drift at 100% may be the most realistic, but more work is clearly required to refine this assumption.
If wave model output is not available then where the waves are in equilibrium with the wind (a 'wind-sea'), it may be a fair approximation to assume that Stokes drift is in the same direction as the wind and proportional to it; e.g. the Doppler sonar observations of Smith (2006) suggest a near-surface Stokes drift of ∼1.25% of the wind speed averaged over the surface 'bubble zone' (∼1.5 m), while the High-Frequency radar observations of Ardhuin et al. (2009) suggest a range 0.5-1.2% of the wind speed for Stokes drift averaged over the top 1 m. However, wave model output is necessary to adequately describe Stokes drift where (1) the wind and wave fields are evolving (2) fetch is limited (near coasts etc), or (3) remotely generated swell is significant (e.g. in the southeastern Pacific; Hanley et al. 2010). Given its importance, Stokes drift should be taken from wave model output where possible.
We also investigated possible numerical errors in the tracking methodologies by carefully undertaking both back-tracking and forward-tracking experiments ( Figure  4). The differences between these two distributions arise from the non-reversibility of the trajectories primarily due to numerical inaccuracies in the tracking scheme in 2-dimensions (which accumulate along the trajectories). These inaccuracies would also apply to other models and other tracking techniques. However, we were here able to examine the net effect of these errors (for the first time in any such study), and found them to be relatively small. In fact, the differences between these figures give a measure of the uncertainty in the method.
Overall, this study shows the utility of our state-ofthe-art assimilative ocean models in providing insights into the drift of particles in the ocean surface. An overview of current and future developments of various assimilative models is described in detail in Schiller et al. (2018). We showed the importance of including Stokes drift, and the need for better quantitative understanding of its magnitude and effect, and also highlighted the high level of accuracy in the numerical schemes which are able to follow the particles over ocean-basin scales (O(1000) km) and gyre timescales (several years). In addition to providing insights into possible aircraft crash sites, the technique also has clear applications for the drift simulation of oil spills, of marine litter and pollution, and of non-swimming larvae, and for studies of the occurrence of invasive species as the ocean currents change in response to climate warming. Note 1. Figure 15 of the CMEMS Quality Information Document, available at http://cmems-resources.cls.fr/ documents/QUID/CMEMS-GLO-QUID-001-024.pdf