Consideration of various aspects in a drift study of MH 370 debris

On March 7, 2014, a Boeing 777-200ER aircraft operated by Malaysian Airlines on the route MH370 from Kuala Lumpur to Beijin abruptly ceased all communications and disappeared with 239 people aboard, leaving a mystery about its fate. The subsequent analysis of so-called satellite 5 ’handshakes’ supplemented by military radar tracking has suggested that the aircraft ended up in the southern Indian Ocean. Eventual recovery of a number of fragments washed ashore in several countries has confirmed its crash. A number of drift studies were undertaken to assist in locating the crash 10 site, mostly focusing either on the spatial distribution of the washed ashore debris or efficacy of the aerial search operation. A recent biochemical analysis of the barnacles attached to the flaperon (the first fragment found in La Réunion) has indicated that their growth likely began in the water of 24°C, 15 then the temperature dropped to 18°C, and then it rose up again to 25°C. An attempt was made in the present study to take into consideration all these aspects. The analysis was conducted by the means of numerical screening of 40 hypothetical locations of the crash site along the so-called 7th 20 arc. Obtained results indicate the likelihood of the crash site to be located between 25.5°S and 30.5°S latitudes, with the segment from 28°S to 30°S being the most promising.


Introduction
On March 7, 2014, approximately 40 minutes after the takeoff, a Boeing-777 aircraft (registration 9M-MRO) operated by Malaysian Airlines (MAS) as MH370 on the route from Kuala Lumpur (Malaysia) to Beijing (China), abruptly ceased all data and voice communications, and disappeared with 239 people aboard, leaving investigators clueless about a possible cause.The subsequent analysis of so-called satellite 'handshakes' (Ashton et al., 2015) supplemented by the primary radar tracking has suggested that the plane turned back, crossed the Malay Peninsula along Malay-Thai border, then flew toward Nicobar Islands in the Strait of Malacca, where it finally turned into the Indian Ocean, as detailed by the Australian Transport Safety Bureau (ATSB, 2014a, b) and the Ministry of Transport Malaysia (MTM, 2017).Eventual recovery of a number of 9M-MRO fragments, which were washed ashore in several countries, has confirmed the crash of the aircraft in the Indian Ocean.A total of 27 suspected and confirmed fragments were found according to the report published by the Malaysian Safety Investigation Team for MH370 (MSIT, 2017) on March 27, 2017: 1 in La Réunion, 2 in Mauritius, 1 in Rodrigues, 6 in Mozambique, 5 in South Africa, 1 in Tanzania, and 11 in Madagascar, with the latest find in January 2017.
Shortly after the disappearance, the Australian Government, whose geographical responsibility for rescue and recovery covers a region of the Indian Ocean where the terminus of 9M-MRO path could have been located according to the satellite data (Ashton et al., 2015), established the Joint Agency Coordination Centre (JAAC, 2014) to assist in the search operation.The Australian Maritime Safety Authority (AMSA, 2014), JAAC and ATSB have conducted an extensive aerial search operation, which lasted from March 18 to April 27, 2014, but failed to locate any debris related to MH370.Although some objects were spotted from the air, subsequent attempts to recover them were unsuccessful.After expiration of the underwater locator beacon, a device emitting acoustic signal to facilitate underwater search, the Australian Government has commissioned an engineering company Fugro N.V., which specializes on marine and geotechnical surveys, to conduct a deep-water highresolution sonar survey of the seabed.The search domain assigned to Fugro N.V. was a band-shaped area in the South Indian Ocean along the so-called 7 th arc from approximately 36°S to 39.5°S latitudes, defined by the ATSB (2015), and then later refined by the Australian Defense Science and Technology Group (Davey et al., 2016).The 7 th arc is a geometric curve on the Earth surface, all points of which are equidistant from the satellite, through which the last 'handshake' was transmitted (Ashton et al., 2015).The actual 7 th arc may slightly differ from the nominal arc due to the uncertainty in the altitude of the aircraft, as well as truncation and measurement errors in the data (e.g., ATSB, 2015; Davey et al., 2016).Despite such an unprecedented effort, the underwater search was unsuccessful, and it was finally called off in January 2017 (MTM, 2017).

10
A number of drift studies were undertaken since then to assist in locating the crash site.Prior to the discovery of the flaperon in La Réunion on July 29, 2015, the studies focused on the analysis of the efficacy of the aerial search, such as García-Garrido et al. (2015).Later, after the flaperon and other 9M-MRO fragments were found, the mainstream approach shifted to the analysis of the probabilities of debris to reach specific locations by known dates starting from various origins along the 7 th arc: the series of numerical and experimental studies undertaken by Griffin et al. (2017) at the Com-20 monwealth Scientific and Industrial Research Organisation (CSIRO); the screening of 25 hypothetical locations along the 7 th arc by Pattiaratchi and Wijeratne (2016) at the University of Western Australia; the numerical modelling conducted by Ormondt and Baart (2015) at Deltares; the study conducted by Maximenko et al. (2015) at the International Pacific Research Center; the study by Jansen et al. (2016) at the Euro-Mediterranean Center on Climate Change; the analysis of drifter trajectories obtained from National Oceanic and Atmospheric Administration's (NOAA) Global Drifter Program (GDP) in relation to MH370 debris by Trinanes et al. (2016).An alternative approach was based on the reverse drift modelling: the studies conducted by the French Government meteorological agency Météo France (Daniel, 2016) and GEOMAR Helmholtz Centre for Ocean Research 35 Kiel (Durgadoo and Biastoch, 2015).The latter, however, did not account for wind forcing, which presumably explains the large difference in its conclusions compared to other studies.
There is an ongoing disagreement between conclusions of these studies with regard to the most likely origin of the debris at the 7 th arc: the latest report published by CSIRO (Griffin et al., 2017) recommends a new search area at around 35°S, backed by the earlier IPRC study (34°S to 37°S).In contrast, Pattiaratchi and Wijeratne (2016) have suggested the crash site to be more likely located between 28.3°S and 33.2°S, narrowing down earlier Jansen et al. (2016) estimates (between 28°and 35°S), being also consistent with Ormondt and Baart (2015).Assuming zero drift angle of the flaperon, Daniel (2016) favors the location north of 25°S, but south of 35°S if the drift angle was set to 18°to the left with regard 50 to wind at the leeway factor of 3.29% -both the parameters experimentally established by the Direction générale de l'Armement (DGA).Griffin et al. (2017) disagrees with these parameters, but confirms observed non-zero drift angles be-tween 0°to 30°, explaining this effect by the longwise asymmetry of the flaperon.
An important feature of the fragment found in La Réunion are barnacles attached to it.Although it was not possible to establish their age, according to De Deckker (2017), who conducted biochemical analysis of the barnacles at the Australian National University, the start of their growth was in the water of approximately 24°C, and then for some time the temperature ranged between 20°C and 18°C, and then it went up again to around 25°C.This additional information has not been previously considered in the drift studies.
Consequently, this study comprises the three major elements to assess the most likely origin of the debris: (1).Efficacy of the aerial search campaign.

Modelling
A total of 40 hypothetical locations of the debris origin along the 7 th arc were screened against the three selection criterion by the means of forward particle tracking technique.The deterministic forcing of each particle in an ensemble was governed by the balance of water and air drag forces, magnitudes of which were assumed to be proportional to the squared relative (with respect to the particle) speeds of the ambient water and air respectively.Surface current velocities were sourced from the Hybrid Coordinate Ocean Model (HYCOM), wind data was sourced from the Global Data Assimilation System (GDAS), as detailed in Section 2.2.2.The stochastic component was modelled using the random walk technique (e.g., Al Rabeh et al., 2000;DHI, 2009;Jansen et al., 2016).Numerical integration was performed in the geocentric Cartesian coordinate system.
All the particles in an ensemble were 'released' from a single starting point at the 7 th arc.Four models with regard to the leeway and drift angle properties of particles were considered.After individual particle tracks were obtained, they were supplemented by respective sea surface temperatures (SST) extracted from the publicly available archives.A subsequent analysis was undertaken in a statistical manner to: (1) estimate the maximum ensemble coverages during the aerial search; (2) assess percentage of particles in ensembles, which could have reached La Réunion by July 29, 2015, and could have been subjected to the temperature variation matching the results of the barnacle analysis; and (3) compare spatial distribution of particles washed ashore with the locations, where MH370 debris were found.

Assumptions
In the frame of this study it was assumed that a particle was subjected to the drag forces induced by water and wind.Parti- cles were assumed to be non-inertial.Impacts of the Coriolis force, Stokes' drift, waves, decay and sinking were neglected (e.g., Kraus, 1972;Al Rabeh et al., 2000).In the local coordinate system (x, y), where the axes x and y correspond to the local west-to-east and south-to-north directions respectively (Figure 1), the location (x i , y i ) of the i-th particle was described by Langevin's equation (DHI, 2009): where (ū i , vi ) is the average (deterministic) velocity of the i-th particle in the local coordinate system, D is the turbulent diffusion term, and ζ, ξ are the random numbers, as detailed in Sections 2.1.3through 2.1.5.

Coordinate systems
On the one hand, velocity of a particle and random walk are formulated in the local coordinate system, where the axis z is normal to the Earth surface.On the other hand, the large study domain dictates the necessity to properly take into consideration the Earth curvature.To perform numerical integration of the governing equations in the geocentric Cartesian coordinate system (X, Y, Z) depicted in Figure 1, where the Earth surface is approximated by WGS'84 ellipsoid with the polar and equatorial axes radii R p =6356752 m and R e = 6378137 m respectively, transformation of coordinates and velocity vectors are required.
The Cartesian coordinates X, Y , and Z of a point on the surface of the ellipsoid described by the longitude ψ and geocentric latitude ϕ can be formulated as: where is the distance between this point and the center of the ellipsoid, as follows from the ellipse equation It should be noted that the backward transformation is required to extract and interpolate surface current and wind data at a location.Respective trigonometric transformations involve solving a 4 th -degree polynomial equation.
The unit vectors − → L , − → M, and − → N, which define directions of 35 the axes x, y, and z of the local coordinate system (Figure 1) are introduced to obtain velocity components in the geocentric system.The outward unit vector − → N = {N X , N Y , N Z , } normal to the surface of the ellipsoid is formulated according to Korn and Korn (1968): where

45
The direction of the axis x is defined by the unit vector − → L : Eq.
(3) and Eq. ( 4) allow for expressing the unit vector − → M = {M X , M Y , M Z } collinear to the axis y in the form of the vector product , so that its components are: 50 Therefore, a velocity vector, the components of which are {ū, v, 0} in the local system, has the following components 55 in the geocentric Cartesian system: Hence, Langevin's equation ( 1) can be formulated in the geocentric Cartesian system as follows: where the relations between the longitude ψ i and geocentric latitude ϕ i of particle's location and its geocentric Cartesian coordinates X i , Y i , Z i is given by Eq. ( 2); the transformation of the velocity components is given by Eq. ( 6); and the local velocity components ūi , and vi are defined from the balance of the deterministic forces.

Deterministic terms
Mathematical description of the dynamics of a floating object is not a trivial problem due to the variety of processes and phenomenon in a near-surface layer, such as surface waves, Stokes drift, flow-object interaction, buoyancy, stratification, etc. (e.g., Kraus, 1972).Such an object is subjected to the dynamic pressure and shear stress forces due to the action of the water and air; its steady-state orientation depends on its buoyancy characteristics and the moments of forces around its principal axes of inertia.Denoting − → ūi = (ū i , vi ) the deterministic velocity of the i-th particle representing a floating object in the local coordinate system, and neglecting the Coriolis force, its motion in the local horizontal plane is described by the equation: where m i is the mass of the object, − − → F w,i and − → F a,i are the average forces caused respectively by water and air flows around the object.The same formulation of these forces as applied in Daniel et al. (2002) and Breivik et al. (2011) to study drift of ship containers, was adopted in this study, although justification for thin nearly-horizontally floating ob-30 jects (Figure 2) could be slightly different.According to the theory of turbulent boundary layer (e.g., Kraus, 1972;Gandin et al., 1955), vertical velocity profiles of the water and air exhibit logarithmic dependence on the distance from the surface , where u * is the friction velocity, κ = 0.41 is von Kármán's constant, and z 0 is the roughness, where the surface is assumed to be non-moving at z = 0.The turbulent shear stresses τ {w,a} = ρ {w,a} u 2 * {w,a} , where ρ {w,a} is the density of the water or air, remains constant through this layer.Hence, the shear stresses acting on the top and bottom surfaces of a thin floating object can be considered proportional to ) are the current and wind velocities at the location of the object at certain reference height and depth (the typical reference height for wind is 10 m above the surface).The dynamic pressures can also be considered proportional to the squared relative velocities of water and air, but at the some other representative distances.Bearing in mind the logarithmic velocity profile, the latter would also be proportional to respectively.Therefore, if the drag forces are collinear to the respective relative velocities of the water and air, the same formulation as in Daniel et al. (2002) would also be applicable for thin objects: where C Dw,i and S w,i are the water drag coefficient and corresponding reference area of the submerged part of the object, C Da,i and S a,i are the air drag coefficient and corresponding reference area of the part exposed to the air, ρ w and ρ a are the water and air densities.Furthermore, Breivik et al. (2011) argued that wave drift forces on small objects (less than 30 m) decay rapidly, and they can be neglected compared to wind forces when the wave length is more than 6 times of the dimensions of a floating object.
It should be noted that composite materials are usually light-weight structures.For example, the recovered flaperon is of approximately 1.6 m x 2.4 m x 0.25 m size and 50 kg weight.Bearing in mind the typical values of the minimum drag coefficient C Da for airfoils are in the range 0.02 to 0.05 with respect to their surface areas, it is easy to show that the drag forces acting on such a fragment of the aircraft would normally be much larger compared to the inertial term in Eq. ( 8), and hence the latter can be approximated by the balance equation: 9), the latter yields: It is easy to see that the solution of Eq. ( 10) is: where the scalar α i (leeway factor) is: For a thin horizontally-floating object, which has equal areas of the surfaces exposed to the seawater and air (S a = S w ), and which is characterized by equal drag coefficients C Da = C Dw , the leeway factor is α ≈3.33%.This theoretical value is in a good agreement with the experimental data for the flaperon (3.29%) estimated by the DGA in the hydrodynamic engineering test facility center in Toulouse (Daniel,10 2016).For a 70%-submerged cube, this formulation yields α = 2.2% assuming the dominance of the dynamic pressure acting on the lateral faces, being also in a good agreement with the experimental factor of 2% (Breivik et al., 2011).
In two of the four considered models (Section 2.2), where the force induced by relative wind was not collinear to its direction, a modified formulation was used instead of Eq. ( 11): where θ is the drift angle, positive counter-clockwise.

Random leeway factor model 20
In contrast to earlier studies (e.g., Daniel, 2016;Pattiaratchi and Wijeratne, 2016;Griffin et al., 2017), an attempt was made in this work to take into consideration the variety of leeway factors, which describe random shapes and flotation characteristics of individual fragments generated by the 25 crash.For this purpose, in one of the four considered types of ensembles (Section 2.2) particles were described by random leeway factors.It was assumed that these particles represented partially submerged thin objects of irregular shapes floating in slightly tilted orientations (Figure 2).For the sake 30 of simplification, it was assumed that the drag coefficients for the air and water are equal (C Da,i = C Ds,i ).Then, according to Eq. ( 12), only the knowledge of the ratio of respective areas k = S a,i /S w,i is required to estimate the leeway factor: In the frame of this study it was assumed that dimensions of individual objects are log-normally distributed.Furthermore, it was assumed that the principal axis of inertia of the object splits it into two parts, areas of which are also log-normally distributed, so that , where {ln S i,1 , ln S i,2 } ∈ N (µ, σ 2 ), and γ ∈ [0, 1] is an independent random parameter to account for the draft of the object (0 -fully submerged; 1 -the center of gravity is at the water surface).Hence the logarithm of the ratio ln(S 1 /S 2 ) = ln S 1 −ln S 2 ∈ N (0, 2σ 2 ) is also normally distributed, with the mean of zero.Here the property was used that the sum of two independent normally distributed random variables is also normally distributed, with its mean being the sum of the means, and its variance being the sum of the two variances (Eisenberg and Sullivan, 2008).Hence, the ratio S 1 /S 2 is log-normally distributed.The modelling was performed assuming ln(S 1 /S 2 ) ∈ N (0, 1).The resulting distribution of the leeway factors of particles in an ensemble in this class of simulations (hereafter referred as the "random leeway" model) is depicted in Figure 3.A random leeway factor assigned to a particle in ensemble was constant in time.

Numerical realization and random walk
The whole integration interval was split into the time steps of ∆t = 15 minutes.Similarly to DHI (2009) particle tracking model, integration of the system of equations (7) for each particle over the time step ∆t was comprised of the two stages, namely deterministic and stochastic: is the velocity of the particle.Unlike DHI (2009) model, which uses the 1-st order discretization method to integrate deterministic terms, the fifth-and sixth-order Runge-Kutta method was used in this work utilizing FORTRAN libraries to ensure particles remain on ellipsoid's surface with sufficient accuracy.During the integration, input current and wind data were bi-linearly interpolated in space, and linearly in time.The vector − → δ i (...) corresponds to a numerical solution for the diffusivity term of the Langevin equation (1).It was treated as a random displacement in the xy-plane locally tangential to ellipsoid's surface.Its components δ x and δ y in the local coordinate system are the random values from the trimmed two-dimensional Gaussian distribution N 2 : where {ζ, ξ} ∈ N 2 (0, 1) are the random numbers, and the standard deviation of the turbulent dispersion σ L = √ 2 D ∆t is assumed to be a function of the horizontal eddy diffusivity coefficient D = D(t, ψ i , ϕ i ).Such a relation between D and σ L was first established by A. Einstein in 1905, who studied diffusion associated with Brownian motion, and since then it was adopted in a variety of random walk models (e.g., DHI, 2009;Jansen et al., 2016).In this study trimming was imposed to discard values δ x and δ y that resulted in displacements exceeding 10 km distance over the time step ∆t.If this criteria was violated, the next pair of random values δ x and δ y was computed.The distance of 10 km was selected as the representative resolution of the ocean circulation model 20 HYCOM, used as a source of the surface current velocities (see Section 2.2.2).
In contrast to Jansen et al. (2016), who applied the constant eddy diffusivity coefficient D = 2 m 2 /s, in this work D was computed according to the well-known Smagorinsky

25
(1963) parameterization, applied in various ocean circulation models (e.g., Blumberg and Mellor, 1987;DHI, 2009): where κ = 0.1 is the constant coefficient (a typical range is from 0.1 to 0.2 according to Blumberg and Mellor (1987)), 30 ∆x and ∆y are the horizontal dimensions of the numerical grid cell, applied for the discretization of the velocity field.A discrete approximation of Eq. ( 16) was used to estimate the eddy diffusivity coefficient D, with subsequent bi-linear interpolation in space and linear interpolation in time at the 35 actual locations of particles.
In the geocentric Cartesian coordinate system the random displacement translates into the 3D vector − → δ i = δX i , δY i , δZ i , components of which are: Such a displacement in the tangential plane causes a particle to move away from the surface of the ellipsoid.However, due to the imposed limitation on the distance, the elevation of the particle does not exceed 8 m after superposition of the random walk.Thus, particle elevations were forced to zeros after applying random walk procedure at each integration time step, while preserving longitudes and latitudes.Box-Muller (1958) transform was used to obtain a pair of pseudo-random numbers {ζ, ξ} ∈ N 2 (0, 1): where τ and ω are the two pseudo-random numbers from the interval (0, 1].To obtain τ and ω, the use was made of the generator developed by Marsaglia and Zaman (1987), claimed to have the period of 2 144 .
Close to a shore, flow velocity was linearly interpolated between zero and velocity in the adjacent cell of the numerical grid of HYCOM.If a particle moved onshore as a result of wind action or random walk, all the subsequent forcing acting on such a particle was nullified, so that it remained at the location where it beached.No specific properties of the shores (e.g.type, slope), or contributing factors such as waves, tides, or storm surges, were considered.
It is worth noting that the transformations between the longitude, latitude and geocentric Cartesian coordinates were performed using extended precision accuracy (80 bits).Conducted tests have shown that the maximum errors in the elevation of a particle arising in a result of a single forwardbackward conversion of the coordinates were of order 0.5 mm using the extended precision compared to 0.5 m using the double precision arithmetic.

Model setup
Four models with respect to the leeway factor and drift angle of particles in an ensemble were considered: (1).Leeway factor of 3.29%; zero drift angle; (2).Leeway factor of 3.29%; drift angle of 18°to the left; (3).Leeway factor of 2.76%; drift angle of 32°to the left; (4).Random leeway factor; zero drift angle.
All the ensembles released at various locations were identical with respect to the properties of particles in these ensembles.Each ensemble comprised 50,000 particles, same as in Pattiaratchi and Wijeratne (2016) study.The second and third models focused specifically on the flaperon path: respective settings corresponded to the flotation characteristics established by the DGA (Daniel, 2016).The fourth model aimed to achieve more realistic representation of the flotation characteristics of the debris generated by the crash by assigning random leeway factors (fixed in time) to the particles of an ensemble; all the 40 ensembles (Section 2.2.1) were identical in terms of the particles they were comprised of.

Screened locations 20
The study domain extended from 20°E to 140°E, 55°S to 15°N.Integration was performed from March 8, 2014 to December 31, 2016 inclusive.The locations of the 40 hypothetical debris origins are summarized in Table 1 and indicated in Figure 4.The coordinates of these locations were estimated from the burst time offset (ATSB, 2015) of the last 'handshake' 00:19, assuming that the aircraft was at the altitude of 10 km.

Model forcing and SST data
The following datasets were used in this study for the model forcing and temperature analysis: • Surface currents were extracted from the Hybrid Coordinate Ocean Model, a data-assimilative isopycnalsigma-pressure coordinate ocean circulation model (Chassignet et al., 2007)  Global Foundation Sea Surface Temperature Analysis, a product developed by the Jet Propulsion Laboratory under NASA MEaSUREs program (JPL, 2015).Spatial resolution: 0.011°x 0.011°; temporal resolution: daily.Detailed information on this data is available at: https://mur.jpl.nasa.gov/and https://podaac.jpl.nasa.gov/dataset/JPL-L4UHfnd-GLOB-MUR.
A direct comparison of the velocity components extracted and interpolated from HYCOM with those of the buoys available from the National Oceanic and Atmospheric Adminis-10 tration's Global Drifter Program (Eliot et al., 2016), which passed through the study domain from March 2014 till December 2016 (a total of 820,801 samples) have shown the RMS errors of 26.6 cm/s and 26.0 cm/s for the easterly and northerly components respectively.However, further analysis is required to understand contribution of wind to these errors (e.g., Griffin et al. (2017) suggests that the average leeway factor of the GDP buoys is around 1.2%), and, more importantly, how they can affect modelling accuracy, particularly with regard to whether they are stochastic or systematic.

Efficacy of the aerial search
An extensive aerial search for MH370 debris, which lasted from March 18 to April 27, 2014, (e.g., JAAC, 2014;AMSA, 2014), failed to find any debris relevant to MH370.Although some suspected objects were observed on March 28 -31 (AMSA, 2014), such as a rectangular object photographed from the specialized Royal New Zealand Airforce (RNAZF) Orion P3K survey aircraft (Figure 5), subsequent attempts to recover them were unsuccessful.Snapshots of the modelled particle locations in the ensembles originated from nine selected locations along the 7 th arc on March 28 and April 5 are depicted in Figure 5 for the random and constant 3.29% leeway factor models.An Animation, which shows computed daily positions of these ensembles and search areas, is presented in the Supplement S1.As seen, the leeway factor had significant influence on the dispersion, which was notably more intense for the ensembles that comprised particles of various leeway factors.Furthermore, difference in the leeway factors could presumably explain the failed attempts to track suspected objects with the help of deployed buoys.It is worth noting that the aerial survey appears to be rather inefficient for the objects characterized by the leeway factor of 3.29% (such as the flaperon), which originated from the locations around 30°S or from the segment from 25.5°S to 28°S of the 7 th arc.
To better understand reasons contributing to the aerial search failure, daily efficacies of the search were analyzed in terms of the percentages E j,k of particles i j,k of each j-th ensemble, which were in the search area Ω k on the k-th day of the search: where N = 50, 000 is the number of particles per ensemble.Respective areas Ω k were obtained by digitizing maps published at AMSA (2014) and JAAC (2014) portals.
As the value of the maximum single-day coverage max k E j,k could be downplayed by other factors affecting detection capability, such as poor weather conditions, the five   6 for the random and constant 3.29% leeway factor models versus origin's latitudes along the 7 th arc.The maximum single-day coverage, respective date of its occurrence, and cumulative coverage over the entire duration of the search campaign are summarized in Table 2 for each screened ensemble, random leeway factor model.The cumulative coverages were defined as the sum k E j,k of daily ensemble coverages over duration of the search campaign.It can be viewed as an additional selection criteria of a more 10 likely origin when peak daily coverages are similar, such as in the case of ensembles No. 21 and 25 (Table 2).As seen, if the crash site was located between 30.5°S and 34.5°S, or between 20°S and 25°S, debris would have been relatively well covered, which would consequently increase chances of their detection.The coverage of MH370 debris would have been notably lower if the crash site was located between 25.5°S and 27.5°S.In particular, for the origins No. 27 and 28, the maximum coverages of approximately 10% occurred more than 6 weeks after the crash, which is a rather 20 low figure bearing in mind that decay and sinking processes were not taken into consideration.A relatively poor coverage is also noted for the origins No. 21 and 23.
Interestingly, a rectangular object spotted from RNZAF Orion P3K on March 28 was in the area, where some debris originating from the location No. 21 (98.07°E,30.03°S) could be expected according to the random leeway model (Figure 5).Furthermore, the timing of the peak coverages (March 28 -31; see Table 2) for the debris originating from this and neighboring locations is consistent with the dates when a number of floating objects were detected there.

SST along the path of the flaperon
One of the major goals of this study is to address a question whether the ambient water temperatures along modelled particle tracks could match the temperature variation derived from the biochemical analysis of the barnacles attached to the flaperon (De Deckker, 2017), and if yes, whether this information could help to further refine the search area.Therefore, those particles were of interest, which satisfied the two following conditions: 1.A particle must have arrived to the Sain-Andre beach, a place where the flaperon was found (assumed coordinates 55.67°E, 20.93°S), by July 29, 2015;  2. The ambient SST at particle's location should have first exceeded 24°C, then dropped down below 18°C, and then risen up again to 25°C or higher.
Due to the inherent uncertainty in the temperature estimations based on the barnacle analysis, a score-based function 5 S was introduced to identify those particles, which approximately satisfied the two aforementioned conditions: Here the term S i,d = exp(−0.5 is responsible for the first condition above; d i = d i (ψ i , ϕ i ) is the ground 10 distance between the i-th particle location on July 29, 2015 and the Sain-Andre beach; d ref = 50 km is the reference distance, which was chosen to be the approximate linear dimension of the Réunion Island.The term S i,θ responsible for the second condition, was formulated in the following way: where θ i (t) is the SST at the i-th particle location at the time t (T s ≤ t ≤ T e , where T s and T e correspond to March 8, 2014 and July 29, 2015 respectively), and The purpose of such a formulation ( 20) is to relax the selection criteria and avoid discontinuities by assigning a positive score to a particle even if it did not arrive precisely to the 10 Sain-Andre beach, or if it did not strictly satisfy the temperature condition (the tolerance allowing for a positive score was set to 1°C).As a result, the maximum score a particle could receive was one.If a particle arrived to La Réunion before July 29, it was still assigned a positive score.If a particle was located at the distance greater than 152 km from the Sain-Andre beach on July 29, it could not receive a score higher than 0.01.If the SST at particle's location never reached 23°C, or never dropped below 19°C after that, such a particle received zero score.

20
The percentages of particles in ensembles, which received scores S>0.01, vs. origin's latitudes along the 7 th arc are shown in Figure 7 for all the four model setups, together with the normal distribution fitting.As seen, for the leeway factors and non-zero drift angles of the flaperon determined by the DGA (Daniel, 2016), the segment centered at 28.2±3°appears to be the most likely area, where the flaperon began its journey.A fraction of particles in the two other models (characterized by zero drift angle), which satisfied the two conditions, could reach a surprisingly high value of 1%, however, both the random leeway factor and constant 3.29% leeway factor models indicated the peak fitted probabilities at 30.8°S and 31.9°Srespectively.The screened origins south of 36.5°S or north of 20°S are deemed to be unlikely starting locations of the flaperon, as none of the four models predicted notable percentage of particles meeting the two requirements above for the corresponding ensembles.
Examples of the first six particle tracks, which received the highest scores S, are depicted in Figure 8 for the scenario corresponding to the leeway factor of 3.29% and drift angle 18°to the left, starting from the origin No. 23 (99.09°E, 28.84°S).All these particles received the scores S i >0.83, and they arrived to La Réunion between June 17 (black track) and July 17 (yellow track), 2015.As seen, the two main reasons for the drop in the ambient water temperature from 23-25°C to as low as 16°C are: 1. Seasonal cooling of the water surface (see comparison of the SSTs on March 8 and June 8, 2014, in Figure 4); 2. Entrapment in counter-clockwise eddies, which could first carry the flaperon northwestward up to 22-25°S latitudes and then southward to 30-33°S (Figure 8).

Beached debris distribution
A total of 27 possibly relevant and confirmed fragments of 9M-MRO were recovered in La Réunion, Mauritius, Rodrigues, Madagascar, Mozambique, South Africa and Tanzania as of April 2017 according to the MSIT (2017); none was ever found in Australia, although a suspected object, the unopened towelette bearing MAS logotype, was discovered at Thirsty Point.Distribution of the found debris offers a useful insight into the possible location of the crash site (e.g., Jansen et al., 2016;Pattiaratchi and Wijeratne, 2016;Griffin et al., 2017), although no consensus was reached up to date with regard to the origin.Similarly to the aforementioned studies, an effort was made in this study to analyze spatial distribution of the washed ashore fragments.
A sample snapshot of the computed particle locations on December 31, 2016 is depicted in Figure 9 for the random leeway model, origin No. 23.Corresponding Animation S2a is included in the Supplement.Considerable beaching of particles of this ensemble was modelled in Africa, Madagascar, La Réunion and Mauritius, but negligible in Australia.The total fractions of particles landed by the ends of 2015 and 2016 are shown in the right top corner for all the screened origins, both the random and constant 3.29% leeway factor models.This result is, in part, due to the West Australian Current, which entrains large percentages of particles from the northern and southern origins (see Animations S3 and S4 in the Supplement), while most of particles from the middle of the screened segment of the 7 th arc remain trapped in the Indian Ocean Gyre.
An interesting result can be obtained by comparing the modelled ratios of the fractions of particles washed ashore in several countries against the ratios of the number of fragments actually recovered in these countries.Such a comparison provides an indication of a number of fragments expected to be found by certain date assuming the same reporting factors.Figure 10 shows computed ratios for the random leeway model, using South Africa as the reference (a total of 5 objects were found there).As seen, more than 5 objects could be expected in Australia, at least 2 objects in Sri Lanka, 9 objects in Tanzania and 3 objects in Kenya, should the origin be south of 35°S.At least one fragment could be expected in Sri Lanka, and several in Kenya and Tanzania for the origins north of 23°S.The best matching segment of the 7 th arc is located approximately between 26.5°S and 31°S.In particular,  the ratio of the percentages of particles beached in Mozambique to those in South Africa was in the range 1.4 to 1.7, being in a good agreement with the actual ratio of 6:5.
To estimate along-shore concentration of MH370 objects expected to be found by certain time, the twodimensional Gaussian smoothing filter was first applied to obtain smoothed concentration of beached particles: where d i is the ground distance between the i-th beached particle and the location of interest (ψ, ϕ), d ref = 5 km is the 10 size of the smoothing filter, and M is the number of beached particles.To compute along-shore density distribution of objects expected to be found, this function was numerically integrated over a relatively narrow band-shaped areas Ω i longwise centered at the shorelines, divided by the respective 15 lengths of the shoreline segments ∆s i , and prorated by the ratio of the number of fragments found in South Africa N SA = 5 to the number of particles landed in South Africa M SA :  An example of the so-estimated concentration of objects is shown in Figure 11  where the unopened MAS towelette was found on July 2 (assumed coordinates 115.3°E, 31°S).A detailed analysis has shown that 8 particles of the respective ensemble landed near Thirsty Point during July 8-11, and 10 more particles arrived before July 28, 2014.All the first 8 particles were characterized by relatively low leeway factors ranging from 0.1% to 0.4%.This suggests a systematic feature rather than a random occurrence.Systematic arrival of particles from the origins north of 27°S or south of 31°S was not predicted earlier than in the last week of July, as seen in Figure 12.The modelling has shown possibility of the sporadic arrival of high-windage particles to the Mossel Bay as early as in February 2015, but the systematic arrival was predicted only in March-May 2016, being consistent with the discovery date.
To understand sensitivity of the results presented in this section to the number of particles in ensembles, a numerical experiment was conducted using 500,000 particles for the origin No. 23, random-leeway factor model (see Animation S2b in the Supplement).The obtained results (Table 3) 10 indicate that the use of 50,000 particles is sufficient to obtain fairly reliable statistics, and further increase of the number of particles would unlikely be beneficial.

Conclusions
The drift study of MH370 debris was conducted by the means of numerical modelling using forward particle tracking technique.A total of 40 hypothetical locations of the crash site along the 7 th arc were screened.The three major aspects were considered: (1) efficacy of the aerial search; (2) ambient water temperatures along the path of the flaperon to La Réunion;

20
(3) the spatial distribution of the washed ashore debris.
The governing equations were numerically integrated in the geocentric Cartesian coordinate system, where the Earth surface was approximated by the WGS'84 ellipsoid.Four models with respect to the leeway factors and drift angles were considered, including a proposed model of random distribution of the leeway factors of particles in an ensemble.
The obtained results indicate significance of the leeway factor in all the three aspects considered.In addition to the uncertainties in the model forcing, assumptions and simplification, judgment about the most likely location of the crash site depends on weights assigned to the aerial search, accuracy of the barnacle biochemical analysis, and probability of fragments not only to be washed ashore, but also recovered and reported.While it does not appear to be possible to confidently point out location of the crash site based on the drift study alone, a few observations can be made with regard to various segments of the 7 th arc: • South of 36°S: Considerable beaching where no fragments were found, particularly in Australia and Sri Lanka; incompatibility with the water temperatures suggested by the barnacle biochemical analysis.
• 34.5 to 36°S: While corresponding areas were poorly surveyed during the aerial search, considerable beaching could be expected in several countries where no fragments were found, particularly in Australia.
• 30.5 to 34.5°S: Excellent aerial coverage of the debris cloud originating from this segment makes the crash site unlikely to be located within it.
• 25.5 to 30.5°S: Consistency with the barnacle temperature analysis; elevated concentration of beached particles, where the fragments of 9M-MRO were found; several 'gaps' in the aerial search; floating objects detected on March 28-31; possible consistency with the early arrival of the MAS towelette to Thirsty Point.
• North of 25.5°S: Inconsistency with the distribution of the washed ashore debris; incompatibility with the barnacle temperature analysis; good aerial coverage of the areas corresponding to the origins from 20°S to 25°S.Summarizing all the above, the most likely area of the crash site appears to be between 25.5°S and 30.5°S, with the segment from 28°S to 30°S being the most promising.This area is consistent with the original definition of a highpriority search zone by the ATSB in June 2014.

Figure 2 .
Figure 2. Schematic representation of a thin floating object.

Figure 3 .
Figure 3. Assumed distribution of leeway factors across an ensemble (constant in time for an individual particle).

Figure 4 .
Figure 4. Locations of the selected hypothetical origins on the 7 th arc, and snapshots of the sea surface temperatures.

Figure 5 .
Figure 5. Snapshots of particle ensembles originating from several locations (indicated with corresponding colors of ensembles) on March 30 for two models, and the areas surveyed till April 2. Sources of the aerial search maps and photograph: JAAC (2014), AMSA (2014).

Figure 6 .
Figure 6.The five largest daily coverages during the aerial search (in the order of descent: black, red, yellow, green and blue bars).

Figure 7 .
Figure 7. Percentages of particles satisfying distance and temperature criterion with the score S>0.01 (LW: leeway factor; DA: drift angle).

Figure 10 .
Figure 10.Expected number of objects to be found in several countries by Dec 31, 2016 vs. origin's latitudes.

Figure 11 .
Figure 11.Expected along-shore concentration of beached MH370 fragments to be found by the end of 2016, origin No. 23.
for the ensemble released from the origin No. 23 (random leeway model); concentrations for all the other studied origins are presented in Figure and Animation S5 in the Supplement.As seen, locations of the elevated concentrations are in a fairly good agreement with the locations, where the fragments were found, except the Rodrigues Island, which was not properly resolved by HYCOM.More fragments could be expected in Tanzania at 7°S and 8.2°S.Screening of the origin No. 23 has revealed that the elevated concentration near Cape Leeuwin in Australia is due to the beaching, which mainly occurred in 2016.During 2014-2015 notable arrival of particles from this origin to Australia took place only around the Windy Harbor and Thirsty Point,

Figure 12
Figure 12 also shows the arrival times of particles to the Mossel Bay in South Africa, where the engine cowling fragment was found.Being the fourth found object, it had covered a long distance over a relatively short time interval.

Figure 12 .
Figure 12.Arrival times of particles to the proximity of the Thirsty Point (Australia) during 2014 and Mossel Bay (South Africa).

Table 2 .
Maximum daily coverage (%), respective date of its occurrence (2014), and cumulative coverage (%) for each screened ensemble, the random leeway model.

Table 3 .
Sensitivity of the percentages of particles washed ashore by Dec 31, 2016 to the number of particles in ensembles (origin No. 23).