A dynamical systems perspective on the absence of debris associated with the disappearance of flight MH 370

Introduction Conclusions References


Introduction
The fate of Malaysia Airlines flight MH370 has been a mystery since its disappearance on the morning of the 8 March 2014. An analysis of radar data, aircraft performance cal- 20 culations and satellite communication (SATCOM) system signaling messages placed the aircraft over an arc, the 7th arc, in the southern part of the Indian Ocean (ATSB, 2014a), where the aircraft's fuel was, presumably, exhausted. Refinements to the analysis of both the flight and satellite data have been continuous since the loss of MH370 and have resulted in the definition of several potential impact regions along the 7th arc, Introduction satellite altimetry, such as AVISO products, which provide mainly geostrophic currents. Similar products integrating surface winds and/or Ekman dynamics are, for instance, OSCAR and SCUD. Alternatively, surface ocean currents can be provided by physics-based computational models which are solved on high-resolution meshes. These provide ocean cir-5 culation currents in space and time on the whole oceanic basin under study and, in order to be representative of the ocean state, they usually assimilate AVISO data and drifters and incorporate other important dynamics and measurements. Examples of these models are the Hybrid Coordinate Ocean Model (HYCOM), or the Operational Mercator global Ocean analysis and forecast system (provided by EU Copernicus). 10 Typically these products provide high resolution space-time data sets but only in areas of interest to the countries funding the consortia providing simulations. In this study the focus is in the middle of the Indian Ocean, where only low resolution data is available.
This study is focused in the analysis of AVISO and HYCOM data predictions for the region and period of interest. Our analysis utilizes fundamental ideas coming from 15 dynamical systems theory applied to Lagrangian transport, which provide sharper insights than the simple drifter tracking supplied by tools such as GNOME or SCUD. The approach taken by the rescue services as described in ATSB (2014a) seems to be closer to a simple drifter tracking approach, and no information is provided on the type of velocity fields used. 20 Dynamical systems theory contributes to this problem by realizing Poincaré's idea of seeking geometrical structures in the phase portrait (for this problem, the ocean surface) that can be used to organize schematically regions corresponding to qualitatively different types of trajectories. These geometrical objects, the stable and unstable manifolds of hyperbolic trajectories, also called Lagrangian structures, are expected to 25 be robust with respect to slight perturbations of the velocity field. The main contribution from this perspective is the consideration that the analysis of individual trajectories may be misleading, and that the right approach to take is to observe the consistency of in-Introduction dividual trajectories with the dynamical barriers highlighted by the Lagrangian skeleton of the velocity field. Finding these structures in the context of geophysical flows is a dynamical systems challenge, as these velocity fields are time dependent and only known on a finite time data set. Their study requires the use of modern techniques in nonlinear dynamical 5 systems theory that achieve a phase portrait representation for systems with these characteristics. In this context many tools have been successful to this end and have provided interesting insights into oceanic problems. For instance, manifolds have been approximated by computing ridges of fields, such as finite size Lyapunov exponents (FSLE) (Aurell et al., 1997) succesfully applied into oceanic contexts (d 'Ovidio et al., 10 2004;Kai et al., 2009) and finite time Lyapunov exponents (FTLE) (Nese, 1989;Shadden et al., 2005Shadden et al., , 2009. Another perspective within the geometrical approach different from Lyapunov exponents is that provided by distinguished hyperbolic trajectories (DHT), a generalization of the concept of fixed point for dynamical systems with a general time dependence (Ide et al., 2002;Ju et al., 2003;Madrid and Mancho, 2009), 15 and their stable and unstable manifolds (Mancho et al., 2006c;Mendoza and Mancho, 2012). In this approach stable and unstable manifolds are directly computed as material surfaces (Mancho et al., 2003(Mancho et al., , 2004. This method has also provided valuable insight into oceanic problems (Mancho et al., 2006a;. Other approaches in this field have been the geodesic theory of Lagragian Coherent Structures 20 (LCS) (Haller and Beron-Vera, 2012), the variational theory of LCS (Farazmand and Haller, 2012), the trajectory complexity measures (Rypina et al., 2011), mesohyperbolicity measures and ergodic partitions (Mezic and Wiggins, 1999;Mezic et al., 2010) and transfer operator methods (Froyland et al., 2012;Froyland and Padberg-Gehle, 2014). 25 A recent tool that reveals phase space structures of general time dependent dynamical systems is the so-called function M, also referred to as Lagrangian Descriptors (see for instance Mancho, 2010, 2012;Mancho et al., 2013), for which several advantages have been discussed in Mancho et al. (2013), Rempel et al. (2013) et al. (2012, 2013) vs. other classical approaches (Shadden et al., 2005;Mezic and Wiggins, 1999). We illustrate the power of this proposed methodology jointly with techniques based on contour advection (Dritschel, 1989;Mancho et al., 2003) to reach conclusions regarding the fate of the plane debris in the weeks after the accident.

5
In this paper we analyze the transport processes arising from AVISO and HYCOM velocity fields. In order to assess the performance (quality) of the ocean state data for the time period and region under consideration we compute for each data set the Lagrangian geometrical skeleton by means of the function M, which is then benchmarked with GDP surface drifter tracks. We find that drifters remarkably follow the Lagrangian skeleton provided by AVISO data. HYCOM data presents high-frequency motions, but still our Lagrangian techniques highlight mesoscale structures analogous to those found in AVISO. By combining contour advection of the presumed impact areas and the geometrical skeleton of the underlying flow we note the presence of ocean mesoscale structures, which could have acted as barriers for the debris. Additionally 15 it is observed that the impact area evolves in such a way that it fills quiet ocean regions, suggesting that they are a very likely destination for the debris. Interestingly, our Lagrangian analysis reaches conclusions that could have guided the search tasks: it highlights regions disregarded by the search efforts, where debris could have been found, and points out that debris could scarcely have visited certain regions subjected 20 to intense search.
The structure of this article is as follows. Section 2 describes the data used in this study. Section 3 describes and discusses the approach to the debris dispersion. Section 4 provides the results and the discussion. Finally, Sect. 5 provides the conclusions.

25
The first data source of velocity fields used for this work has been obtained from satellite altimetry produced by Ssalto/Duacs and distributed by AVISO, with support from Introduction CNES, and are downloadable from http://www.aviso.oceanobs.com/duacs/. In particular we have considered the product UPD, which is reconstructed from four satellites. The data are given on a spatial grid of 1080 × 915 points (longitude/latitude). The latitude ranges from 82 • S to 81.9746 • N using a Mercator projection, and the longitude ranges from 0 to 359.667 • using a uniform grid. The spatial precision is thus 1/3 • at 5 the Equator and is provided with daily frequency. AVISO data essentially reproduces geostrophic currents. A priori, the convenience of using this type data for our study is supported by work discussed in Olascoaga et al. (2013) and Mendoza et al. (2014), and eventually is confirmed by the results reported in the discussion section. A second source used for the velocity fields is the Hybrid Coordinate Ocean 10 Model (HYCOM) (Bleck, 2002;Chassignet et al., 2007), which provides velocity fields throughout the entire ocean depth. This model assimilates satellite sea-surface height and temperature data by means of the Navy Coupled Ocean Data Assimilation (NCODA) system and it is forced by surface winds and air-sea fluxes. In particular, we have considered the GLBa0.08 data set (experiments expt_91.0 and expt_91.1), 15 which can be downloaded from http://hycom.org/dataserver/glb-analysis/. These experiments use 32 layers along the vertical coordinate and are carried out over a global Mercator curvilinear grid for the ocean, where latitude ranges from 78 • S to 47 • N and the equatorial resolution is 1/12 • , resulting in a grid point spacing of approximately 7 km on average. Regions north of 47 • N are represented with a bipolar patch. The 20 data is provided daily. HYCOM data incorporates several physical factors that introduce high-frequency variability, and produce velocity fields that, at least in the very upper layers, are far from geostrophy. In the next section we will discuss issues of finite resolution and noise effects in relation to the debris dispersion problem. In the calculations performed we 25 attempt to filter the high-frequency variability by using the HYCOM data layer at 50 m depth. A similar depth choice for this model is used, for instance, in Sulman et al. (2013) or for CUPOM (Colorado University Princeton Ocean Model) in Branicki and Kirwan (2010). In the latter work the authors show how the resulting Lagrangian struc- ture is close to 2-D surfaces that extend nearly vertically into the water column. The full 3-D Lagrangian structures can be obtained by a vertical extension of the evolving stable and unstable manifolds calculated in the 2-D plane, but this is valid only throughout a layer whose thickness depends on the ratio of the characteristic velocity within the 2-D slice to the characteristic vertical shear of the horizontal velocities (Branicki and 5 Kirwan, 2010; Branicki et al., 2011). We confirm that this is the case in our study. Finally, in order to benchmark the Lagrangian structures calculated for the velocity field data sets described above, surface drifter tracks are obtained from data distributed by the Global Drifters Program (GDP), from the National Oceanic and Atmospheric Administration (NOAA) and Atlantic Oceanographic and Meteorological Laboratory (AOML). These tracks can be downloaded from http://www.aoml.noaa.gov/phod/ dac/index.php. Each drifter has an identification number with five digits, which allows for the identification of different drifters. The sampling of the drifter position is once every 6 h. 15 The debris produced after a plane's accident is due to the breakup of the plane and it depends on how it enters the water. Large heavy and unbroken pieces of fuselage would sink rapidly. However, there exist reports of plane accidents (Chen et al., 2015) which produced debris spread over a wide area, with light pieces that might have floated for a long time. This article, as described in the Introduction, is focused on 20 the search for this type of debris, such as flat wide objects not surpassing the waterline, and thus not subjected to wind sailing effects, that would have been mainly driven by ocean currents. We will assume a similarity between this debris and drifter motion. Kuznetsov et al. (2002) have examined the role played by Lagrangian structures obtained in the Gulf of Mexico. In particular, they have analyzed data from the Colorado Recent studies (Tallapragada and Ross, 2008;Beron-Vera et al., 2015) have considered deviations in the evolution of drifters from that of purely advected particles by taking into account inertial effects produced by the buoyancy of the object and their finite size. Their approach considers the Maxey-Riley equation (Maxey and Riley, 1983) which holds for small rigid spheres. Considering the diversity of potential shapes for 10 the debris objects, and that mainly we are considering flat, wide objects, quite different from a sphere, we do not consider this approach in our study. We disregard this choice because it is at least as speculative as the straightforward Lagrangian approach for purely advected particles, and for this reason the latter is the one we take. In the end, the agreement found between the GDP drifters tracks and the mesoscale features that 15 our tools highlight supports this decision. In this way we consider the motion of the airplane debris as that of purely advected particles. Such particles follow trajectories x(t) in the ocean that evolve according to the dynamical system:

Modelling debris dispersion
where v(x(t),t) is the velocity field of the ocean region of interest. In our analysis, as 20 previously explained, we will assume that the motion of particles is mainly horizontal. The equations of motion that describe the horizontal evolution of particle trajectories on a sphere of radius R are: Here the variables (λ, φ) are longitude and latitude, u and v respectively represent the eastward and northward components of the velocity field provided by the data set. For our Lagrangian analysis, which is based on Lagrangian Descriptors (function M) and also on contour evolution, it is necessary to integrate Eqs. (2) and (3) in order to compute particle trajectories. As the velocity field for the system (Eqs. 2 and 3) 5 is provided solely on a discrete space-time grid, the first issue to deal with is that of interpolation. We have used bicubic interpolation in space and third order Lagrange polynomials in time according to the details given in Mancho et al. (2006b) and Mendoza et al. (2014). This assumes smoothness of the velocity field below the resolution which, on the other hand, is rather low. 10 An open question here is how robust are Lagrangian techniques when dealing with ocean motions that do not vary smoothly at small scales, as those having an inherent turbulent nature. Recent work by Hernández-Carrasco et al. (2011) analyze this question by looking at the impact on FTLE of the two major effects when dealing with real data, namely, noise and dynamics of unresolved scales. Their work concludes that, 15 although some features are lost, this tool still gives an accurate picture of the oceanic transport properties. We take these conclusions as supportive of our assumptions in the interpolation scheme that we use. On the other hand, effects of high frequency variations on the velocity field are also included in the HYCOM data. In the last section we will highlight common features between this and the more smooth AVISO data, which 20 will support also the conclusions in Hernández-Carrasco et al. (2011).

The Lagrangian skeleton
The global nature of particle trajectories generated by Eq. (1) Madrid and Mancho, 2009;Mancho et al., 2013), defined as follows: Here || || stands for the modulus of the velocity vector. At a given time t 0 , function M(x 0 , t 0 ) measures the arclength traced by the trajetory starting at x 0 = x(t 0 ) as it evolves forwards and backwards in time for a time period τ. The structure of the function M shows, at low τ values, a smooth pattern such as that visible in Fig. 1a. There, M 10 has been evaluated for AVISO data, on the date of the plane crash and near one of the search areas, using an integration period of τ = 5 days. On the other hand, Fig. 1b (computed for τ = 20 days) illustrates how the structure of M evolves for large τ towards less regular structures. By this we mean that sharp changes of M values occur in narrow gaps, forming filaments that highlight stable and unstable manifolds. A thorough 15 explanation of this effect is discussed in Mancho, 2010, 2012), Mancho et al. (2013) and Lopesino et al. (2015). Figure 1c shows a different projection for the latter case, which highlights the singular features of M aligned with stable and unstable invariant manifolds, and their crossing in a hyperbolic point. We remark that the information contained in these figures is obtained by means of integration of particle 20 trajectories. In fact some of the figures reported in this article use tens of thousands of them. The difference in this analysis with respect to that performed by the drifter tracking approach taken by ATSB (2014a, b), or other tools such as GNOME or SCUD, is the way in which this information is displayed. The latter tend to represent this output as spaghetti diagrams, while our Lagrangian perspective finds geometrical structures Large M values, represented in reddish colors, are related to regions of high speed fluid (such as jets), while bluish colors denote calm regions. Figure 1d schematically shows the manifold skeleton and how particles, such as drifters or plane debris would 5 evolve in the neighborhood of this structure. Invariant manifolds are dynamic barriers that cannot be crossed by purely advected particles and control transport in the region. The cross point of the stable and unstable manifolds, the hyperbolic point, drastically bends the trajectories. Tracer blobs evolve constrained by the manifold structure. Figure 1e illustrates this point by showing elongation of blob contours along the unstable 10 manifold of a hyperbolic trajectory and squeezing along its stable manifold. It is important to remark here that, from the dynamical system perspective, the main distortion agents for the contours are those causing their filamentous shape and, as explained here, these agents are the hyperbolic trajectories and their stable and unstable manifolds. Additionally, in the flow there exists regions dominated by trajectories with elliptic 15 type stability, in which, contrary to what happens near hyperbolic trajectories, blobs remain coherent instead of forming filaments. This is illustrated in Fig. 1f. Oceanic eddies are an example of regions displaying such a behaviour (see, e.g. Chelton et al., 2011) and this is also the case for oceanic jets, in which blobs are transported at a high speed but not distorted into filaments (see e.g. Wiggins and Mancho, 2014). Therefore, the in-20 terplay between stretching and confinement is an essential ingredient of fluid transport processes.
In this work, the importance of how water masses evolve is clear because of the uncertainty in the plane's impact point, and thus, in order to reach conclusions on the fate of the plane debris, it is important to track the time evolution of areas in which the plane 25 potentially could have crashed, and recognize the constraints on these regions found in the ocean during the period after the accident. The evolution of these areas is tracked with the contour advection algorithm developed by Dritschel (1989), but including some modifications explained in Mancho et al. (2003Mancho et al. ( , 2004Mancho et al. ( , 2006c

Analysis of the MH370 plane debris dispersion
The goal of this section is to analyze the search strategy followed by the Australian Maritime Safety Authority (AMSA), in order to see if, in the light of the tools and data described in the previous sections, alternative search strategies could have been suggested. 5 One of the major challenges in the first stages of the search was to know how debris had scattered in an always moving ocean. On Tuesday 18 March and on Wednesday 19 March search areas were defined by AMSA according to possible flight paths deduced from satellite and aircraft performance data which assumed impact points along the 7th arc. The charts released by AMSA for the different surface search areas 10 are available at http://www.amsa.gov.au/media/incidents/mh370-search.asp. Figure 2a shows the potential flight paths in white, and the 7th arc in black. The presumable impact area along the 7th arc is depicted in pink. Its width, which is at around 93 km, expresses uncertainties in the impact point. Our discussion is based on this uncertainty band, which was considered as the maximum priority in the reports by ATSB 15 (2014a, b). However, we note that the satellite analysis also has highlighted a wider search band around this narrower strip.
The search area released by AMSA during these days appears in the figure surrounded by a dashed line. Moreover, the gray shaded contour shown in this figure represents the advection of the presumed initial impact area until the 19 March, according 20 to ocean velocity data distributed by AVISO. This analysis thus suggests a smaller search area than the one determined by AMSA for the search operations that took place during the 18 and 19 March. However, one remaining question that needs to be addressed in order to support this conclusion is, how accurate are AVISO data when representing the true ocean state? We provide some discussion of this point in terms of 25 the tools described in the previous subsection. We observe that the advected contour stretches along unstable manifolds, which act as the main deforming agents for this blob. Stable and unstable manifolds and eddies are recognized in the background of  Figure 3a and b proceeds with a similar analysis to that of Fig. 2a and b, but using HYCOM data instead. In the HYCOM simulations, layers near the surface are strongly dominated by wind forcing and, as a result, the velocity fields are subjected to high frequency variations. As explained in Sects. 1 and 2, our study assumes a similarity between debris motion and that of drifters. Mancho et al. (2006a), Branicki and Kirwan 25 (2010) and Sulman et al. (2013) justify the study of Lagrangian structures from velocity fields by considering layers just below the very upper layers which are slightly wind filtered. In Fig. 3a and b the background coloring shows 2-D Lagrangian structures obtained at 50 m depth. The Lagrangian structures now are much more blurred than Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | those obtained from AVISO, but still a comparison is feasible in terms of the mesoscale structures highlighted by the M function. HYCOM shows mesoscale features linked to the circumpolar current further spread into the north, and the calm region close to drifter 16545 still seems to be present, but rather diffuse and displaced south of the drifter. Figure 3a illustrates that the spreading of the presumed impact area is broader 5 than that computed from AVISO data, but still smaller than the AMSA search area. Further details are found in movie S2. The calm area in which drifter 16539 seemed to be anchored, according to Fig. 3b, now appears as a more active region with diffuse structures. However the HYCOM analysis likewise confirms that the debris of a plane crashing along the 7th arc could have barely reached this region in this period of time. 10 On the 28 March further refinements of the flight path shifted the search area up to the north, and the newly determined most likely impact area was then the one depicted in Fig. 4a ATSB (2014a) in the northeast direction above the Broken Ridge along the 7th arc. From the 28 March until the 3 April the search services were focused on the area surrounded by the dashed line in Fig. 4a. The background coloring of this figure   15 shows the Lagrangian structures obtained from the M function evaluated from AVISO data on the 3 April 2014. The presumed initial impact area is advected until this day and represented in the gray shaded region. Its filamented structures clearly tend to be aligned along the unstable manifolds visible in the singular features of M. Numerous drifters are represented with magenta dots and they move consistently with the high-20 lighted Lagrangian skeleton. This is further confirmed in the video S3. For instance drifter 56566, which, on the day of the accident, is within a strong mesoscale structure, an oceanic eddy placed at the east of the presumed impact area, remains within the eddy during several months, being transported with it towards the west. Similarly, drifters 56947 and 56512 are trapped at some point in weaker eddies. Drifters 56557 The strong eddy containing the drifter 56566, behaves like a capsule where surrounding waters cannot penetrate. The advected gray shaded region deforms around it, showing it is a forbidden region for the debris. On the 3 April, this large vortex invaded the search area, and thus our analysis suggests that the region occupied by this moving coherent structure could have been removed from the region where the search 5 efforts were focused.

NPGD
Another interesting mesoscale element in this area is the strong jet which crosses the presumed impact area at the north end. The feature is visible in orange-reddish colors and its presence is confirmed by the motion of drifter 56554 which navigates towards the east between the 27 March until the 11 April. The jet persists since the 10 impact day (the 8 March) until at least the end of May and as we will discuss below is a dynamic barrier. It separates the blob located in the impact area into two pieces which remain unmixed for the whole period of study, a small one at its north side and the larger one at south side.
Between the 3 and the 11 April the search area was moved up to the north, to regions 15 rather diverted from the probable impact area displayed in Fig. 4a. We ignore this part and we pursue our discussion with the search targets held from the 12 to 28 April, which are consistent with the impact area under consideration. The search areas after the 12 April are reported by the Australian Government (see at http://www.jacc.gov.au/ media/releases/2014/april/index.aspx), and correspond to the regions surrounded by 20 dashed lines in Fig. 4b. There, it is shown again the presumed impact area along the 7th arc, and the advection of this contour by the 15 April. From the background coloring of the M function, we conclude that as time evolves the contour tends to form filaments aligned with the unstable manifolds and accumulates into calm regions which gather most of the material (see also the video S3). From the figure it is clear that most of 25 the impact area has scattered into the calm region below the strong jet described in the previous paragraph. Interestingly, this region was left out of the search area for this period, and reversely the search area barely contains material coming from the impact area. As time evolves the gray shaded contour evolves and penetrates more and more into the search panel at the right, but just at the south of the strong jet. We see that the jet acts as a barrier which prevents most of the material of the likely impact area from being transported towards the north, this had the effect of paradoxically concentrating the main search in an area where, according to our analysis, there was likely no debris. Figure 5a and b shows a similar analysis to that of Fig. 4a and b using HYCOM data 5 at 50 m depth. Mesoscale structures in this data are consistent with those of AVISO. For instance, the large eddy which isolates drifter 56566 is present in a similar position, invading the search area from the 28 March to the 3 April, confirming again that this part could have been excluded from the search. The jet structure splitting the impact region in two is also present at the same position, although somewhat perturbed by 10 wind effects. Similarly to AVISO, much of the material contained in the impact area would evolve filling calm regions which, roughly speaking, are positioned the same. The analysis of Fig. 5b is similar to that of Fig. 4b, and it says that, according to our diagnosis, on the 15 April most of the material in the impact area is outside the search area. Video S4 completes these conclusions.
By the 8 October 2014 the Australia Transport Safety Bureau ATSB (2014b) released a flight path analysis update which indicated that the next underwater phase of the search should be prioritized further south along the 7th arc, and thus a detailed bathymetry was required, and efforts were initiated to obtain this for the zone depicted in Fig. 6. A Lagrangian analysis of the horizontal motions at 2000 m depth confirms that 20 at these depths, the currents are more calm than in upper layers, and consequently the evolution of the impact contour does not spread so much at depth in the ocean. On the other hand, this phase of the search is intended for heavier objects, whose motion would be more gravity dominated than current dominated, and thus horizontal dispersion should not be large for these objects. Movie S5 and Fig. 6 confirm that the north 25 and south bounds of the underwater impact area are the regions in which dispersion is the largest.

Conclusions
This work addresses, from the dynamical systems perspective, the fate of the MH370 plane debris after its presumed crash on the morning of the 8 March 2014. This point of view provides a powerful tool for dealing with the difficulty arising from the uncertainty of predictions, which comes from the fact that there are several data sets that de-5 fine multiple ocean states. The findings described in this paper highlight several facts that could have been relevant to the search for debris associated with MH370. First, we have shown that, according to several available data sources, impact debris from the MH370 plane would have been mainly scattered into calm ocean regions since the material from the likely impact areas tends to fill these areas. Second, mesoscale 10 structures such as ocean eddies and jets may have played a key role in determining forbidden regions and barriers for debris. The consistency found between data obtained from different and independent sources (AVISO, HYCOM, drifters) back the evidence that the mesoscale structures we have described were present in the Indian Ocean during the accident time frame, and that considering these structures could have guided 15 an efficient search procedure. These facts also back the reliability of our dispersion analysis, which would have suggested channeling efforts to regions disregarded from scheduled search areas and bypassing certain regions that were subjected to intense search. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | de la Cámara, A., Mancho, A. M., Ide, K., Serrano, E., and Mechoso, C.: Routes of transport across the Antarctic polar vortex in the southern spring, J. Atmos. Sci., 69, 753-767, 2012. 1202 de la Cámara, A., Mechoso, R., Mancho, A. M., Serrano, E., and  Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Ide, K., Small, D., and Wiggins, S.: Distinguished hyperbolic trajectories in time-dependent fluid flows: analytical and computational approach for velocity fields defined as data sets, Nonlin. Processes Geophys., 9, 237-263, doi:10.5194/npg-9-237-2002, 2002