Message Boards Message Boards

Ocean currents: from Fukushima and rubbish, to Malaysian airplane MH370

In this post I will use data from NASA's ECCO2 project (http://ecco2.jpl.nasa.gov/) to simulate various scenarios: the movement of radioactive particles from the Fukushima nuclear power plant and the accumulation of rubbish in the oceans. I download about 20 years worth of oceanographic data and study how the flow of water might transport particles. A slight modification of the program might be useful to determine the crash site of Malaysian airplane MH370. The work was done with Bjoern Schelter, who is member of this community.

We will first download the vector fields for u and v direction. We obtain data for different depth, but I will only use surface currents. We first generate a list of filenames.

filenamesU = ("http://ecco2.jpl.nasa.gov/opendap/data1/cube/cube92/lat_lon/quart_90S_90N/UVEL.nc/UVEL.1440x720x50." <> # <> ".nc.ascii?UVEL[0:1:0][0:1:0][0:1:719][0:1:1439]" & /@ (StringJoin[{ToString[#[[1]]], StringTake[StringJoin[ToString /@ PadLeft[{#[[2]]}, 2]], -2], StringTake[StringJoin[ToString /@ PadLeft[{#[[3]]}, 2]], -2]}] & /@ Table[Normal[DatePlus[DateObject[{1992, 1, 2}], 3*n]], {n, 0, 2556}]));

filenamesV = ("http://ecco2.jpl.nasa.gov/opendap/data1/cube/cube92/lat_lon/quart_90S_90N/VVEL.nc/VVEL.1440x720x50." <> # <> ".nc.ascii?VVEL[0:1:0][0:1:0][0:1:719][0:1:1439]" & /@ (StringJoin[{ToString[#[[1]]], StringTake[StringJoin[ToString /@ PadLeft[{#[[2]]}, 2]], -2],StringTake[StringJoin[ToString /@ PadLeft[{#[[3]]}, 2]], -2]}] & /@ Table[Normal[DatePlus[DateObject[{1992, 1, 2}], 3*n]], {n, 0, 2556}]));

These filenames cover a time interval from the beginning of 1992 to 2012. We then download the actual data (which can take really long):

Monitor[For[k = 1, k <= Length[filenamesV], k++, 
  Export["~/Desktop/OceanVelocities/velV" <> ToString[k] <> ".csv", 
   ToExpression[StringSplit[StringSplit[#, "],"][[2]], ","]] & /@ 
    StringSplit[Import[filenamesV[[k]]], "VVEL.VVEL[VVEL.TIME="][[
     2 ;;]] ]], 
 ProgressIndicator[Dynamic[k], {0, Length[filenamesV]}]]

and

Monitor[For[k = 1, k <= Length[filenamesU], k++, 
  Export["~/Desktop/OceanVelocities/velU" <> ToString[k] <> ".csv", 
   ToExpression[StringSplit[StringSplit[#, "],"][[2]], ","]] & /@ 
    StringSplit[Import[filenamesU[[k]]], "UVEL.UVEL[UVEL.TIME="][[
     2 ;;]] ]], 
 ProgressIndicator[Dynamic[k], {0, Length[filenamesU]}]]

This will download about 16GB worth of data. You will, however, get reasonable results if you only use 30 frames, i.e. if you substitute Length[filenamesU] and Length[filenamesV] by 30. We will now create a figure for the first frame.

velV1 = Import["~/Desktop/OceanVelocities/velV1.csv"];
velU1 = Import["~/Desktop/OceanVelocities/velU1.csv"];

We can now generate the velocity vectors for each point on the surface of the earth

veltot = Table[{velU1[[i, j]], velV1[[i, j]]}, {i, 1, 720}, {j, 1, 1440}];

and then clean the data

veltot2 = Table[If[((StringQ[#[[1]]] || StringQ[#[[2]]]) & /@ {veltot[[i, j]] })[[1]], {0., 0.}, {velU1[[i, j]], velV1[[i, j]]}], {i, 1, 720}, {j, 1, 1440}];

We can also calculate the speed at each point.

veltottot = Table[If[((StringQ[#[[1]]] || StringQ[#[[2]]]) & /@ {veltot[[i, j]] })[[1]], 0., Norm[{velU1[[i, j]], velV1[[i, j]]}]], {i, 1, 720}, {j, 1, 1440}];

This next function will produce a visualisation of the speed profile for the first frame.

Show[ArrayPlot[Reverse[veltottot], PlotLegends -> Automatic, ImageSize -> Full, ColorFunction -> "TemperatureMap"], ListVectorPlot[Transpose[veltot2], VectorPoints -> 50]]

enter image description here

Note the incredible fine-structure in that image. We could now produce an animation for all frames. This will, however, be too large for the memory of most computers, so I would not execute it. If you only use the first 30 frames you get a reasonable idea of its working.

framesflow = {}; Monitor[
 For[k = 1, k <= Length[filenamesU], k++, 
  velU1 = Import["~/Desktop/velU" <> ToString[k] <> ".csv"]; 
  velV1 = Import["~/Desktop/velV" <> ToString[k] <> ".csv"];
  veltot = Table[{velU1[[i, j]], velV1[[i, j]]}, {i, 1, 720}, {j, 1, 1440}]; 
  AppendTo[framesflow, ArrayPlot[Reverse[Table[If[((StringQ[#[[1]]] || StringQ[#[[2]]]) & /@ {veltot[[i,j]] })[[1]], 0., Norm[{velU1[[i, j]], velV1[[i, j]]}]], {i, 1, 720}, {j, 1, 1440}]], PlotLegends -> Automatic, ImageSize -> Full, ColorFunction -> "TemperatureMap"]]], k];

This following command exports all the frames as individual gif files, which can later be combined into a gif animation.

Monitor[For[k = 1, k <= Length[filenamesU], k++, 
   velU1 = Import[
     "~/Desktop/OceanVelocities/velU" <> ToString[k] <> ".csv"]; 
   velV1 = Import[
     "~/Desktop/OceanVelocities/velV" <> ToString[k] <> ".csv"];
   veltot = 
    Table[{velU1[[i, j]], velV1[[i, j]]}, {i, 1, 720}, {j, 1, 1440}]; 
   Export["~/Desktop/MovieOut/frame" <> ToString[1000 + k] <> ".gif", 
    ArrayPlot[Reverse[Table[If[((StringQ[#[[1]]] || StringQ[#[[2]]]) & /@ {veltot[[i,j]] })[[1]], 0., Norm[{velU1[[i, j]], velV1[[i, j]]}]], {i, 1, 720}, {j, 1, 1440}]], PlotLegends -> None, ImageSize -> {1440, 730}, ColorFunction -> "TemperatureMap"]]], k];

Another advantage of saving the frames is that we can use them as background for various scenarios and thereby decrease our computation time substantially. Here is a short animation of the first couple of frames.

enter image description here

Now we get the trajectories of radioactive particles that are released into the ocean at the approximate site of the Fukushima power plant. We make the assumption that the particles simply follow the flow. Also, we use the vectorfield starting on 2nd January 1992. The overall patterns of the ocean currents seem to be reasonably stable over the years so that for this conceptual study this assumption will have to suffice. It is of course easy to adapt to the actual date.

We first need to introduce a formula to convert the data grid to the surface of the earth.

ClearAll["Global`*"];
f[teta_, r_, vx_, vy_] := {vx/(r*Sin[2*Pi*teta/4./360.]*2*Pi/360/4)*24*3600, vy/(N[r*2*Pi/360/4])*24*3600};
teta = 1; r = 6360000;

We now place 20000 particle into the see close to the location of Fukushima. We add a small random number to the position to model a "cloud" of particles.

teilreinpos = 
 Table[{3, 511 + RandomVariate[NormalDistribution[0., 1.]], 568 + RandomVariate[NormalDistribution[0., 1.]]}, {m, 1, 20000}]; 

This is the main part of the program. We use the velocity to update the position of the particles. We then use the back-ground frames generated above and overlay the particle positions.

trajectories = {teilreinpos};
Monitor[Do[
   velVt = Table[
     If[StringQ[#], 
          0.1*(ToExpression[#[[
                 3]]]*10^(ToExpression[#[[1]]]) ) & /@ { 
            StringSplit[#, {" ", "*"}]}, 1.*#] & /@ # & /@ 
      Import["~/Desktop/OceanVelocities/velV" <> ToString[k] <> 
        ".csv", "CSV"], {k, i - 2, i + 2}];
   velUt = 
    Table[If[StringQ[#], 
          0.1*(ToExpression[#[[
                 3]]]*10^(ToExpression[#[[1]]]) ) & /@ { 
            StringSplit[#, {" ", "*"}]}, 1.*#] & /@ # & /@ 
      Import["~/Desktop/OceanVelocities/velU" <> ToString[k] <> 
        ".csv", "CSV"], {k, i - 2, i + 2}];
   xfunc = ListInterpolation[velUt, {{1, 5}, {0, 720}, {0, 1440}}];
   yfunc = ListInterpolation[velVt, {{1, 5}, {0, 720}, {0, 1440}}];
   AppendTo[trajectories, 
    Nest[({#[[1]], Mod[#[[2]], 720], 
           Mod[#[[3]], 1440]}) & /@ (# + {0.1, 
            0.1*f[#[[2]], r, xfunc[#[[1]] - i + 3., #[[2]], #[[3]]], 
               yfunc[#[[1]] - i + 3., #[[2]], #[[3]]]][[2]] , 
            0.1*f[#[[2]], r, xfunc[#[[1]] - i + 3., #[[2]], #[[3]]], 
               yfunc[#[[1]] - i + 3., #[[2]], #[[3]]]][[1]] } &) , # ,
        10] & /@ trajectories[[-1]]]; 
   Export["~/Desktop/OilSpillFrames/frame" <> ToString[1000 + i] <> 
     ".gif", Show[
     Import["~/Desktop/MovieOut/frame" <> ToString[1000 + i] <> 
       ".gif"], 
     Graphics[{Red, Disk[#, 1]} & /@ (0.9361111 # + {44., 28.} & /@ 
         trajectories[[i - 2, All, {3, 2}]]), 
      PlotRange -> {{0, 1440}, {0, 720}}]]];
   , {i, 3, 1238}], i];

Here are a couple of frames:

enter image description here

All frames yield the following animation.

enter image description here

Also see the higher quality animation on Youtube. https://youtu.be/8icgMg6lt0Y

Similarly we can study where rubbish, such as plastic bags would accumulate in the oceans due to the currents; these are also called gyres of marine debris particles. See also this link.

We start out just like for the Fukushima simulation:

ClearAll["Global`*"];
f[teta_, r_, vx_, vy_] := {vx/(r*Sin[2*Pi*teta/4./360.]*2*Pi/360/4)*24*3600, vy/(N[r*2*Pi/360/4])*24*3600};
teta = 1; r = 6360000;

In this case we have to distribute the particles randomly in all oceans. So I will basically distribute them everywhere and then ingore the ones over land mass, by checking whether the background flow speed is zero. So we first calculate the speed.

velVt = Table[
   If[StringQ[#], 
        0.1*(ToExpression[#[[3]]]*10^(ToExpression[#[[1]]]) ) & /@ { 
          StringSplit[#, {" ", "*"}]}, 1.*#] & /@ # & /@ 
    Import["~/Desktop/OceanVelocities/velV" <> ToString[k] <> ".csv", 
     "CSV"], {k, 3 - 2, 3 + 2}];
velUt = Table[
   If[StringQ[#], 
        0.1*(ToExpression[#[[3]]]*10^(ToExpression[#[[1]]]) ) & /@ { 
          StringSplit[#, {" ", "*"}]}, 1.*#] & /@ # & /@ 
    Import["~/Desktop/OceanVelocities/velU" <> ToString[k] <> ".csv", 
     "CSV"], {k, 3 - 2, 3 + 2}];

xfunc = ListInterpolation[velUt, {{1, 5}, {0, 720}, {0, 1440}}];
yfunc = ListInterpolation[velVt, {{1, 5}, {0, 720}, {0, 1440}}];

We then distribute the particles everywhere

teilreinpospre = Table[{3, RandomReal[{0, 720}], RandomReal[{0, 1440}]}, {k, 1, 30000}];

and delete the particles with zero-speed-background.

teilreinpos = Select[teilreinpospre, xfunc[3, #[[2]], #[[3]]] != 0. && yfunc[3, #[[2]], #[[3]]] != 0. &];

Now we iterate as before.

trajectories = {teilreinpos};
Monitor[Do[
   velVt = Table[
     If[StringQ[#], 
          0.1*(ToExpression[#[[
                 3]]]*10^(ToExpression[#[[1]]]) ) & /@ { 
            StringSplit[#, {" ", "*"}]}, 1.*#] & /@ # & /@ 
      Import["~/Desktop/OceanVelocities/velV" <> ToString[k] <> 
        ".csv", "CSV"], {k, i - 2, i + 2}];
   velUt = 
    Table[If[StringQ[#], 
          0.1*(ToExpression[#[[
                 3]]]*10^(ToExpression[#[[1]]]) ) & /@ { 
            StringSplit[#, {" ", "*"}]}, 1.*#] & /@ # & /@ 
      Import["~/Desktop/OceanVelocities/velU" <> ToString[k] <> 
        ".csv", "CSV"], {k, i - 2, i + 2}];
   xfunc = ListInterpolation[velUt, {{1, 5}, {0, 720}, {0, 1440}}];
   yfunc = ListInterpolation[velVt, {{1, 5}, {0, 720}, {0, 1440}}];
   AppendTo[trajectories, 
    Nest[({#[[1]], Mod[#[[2]], 720], 
           Mod[#[[3]], 1440]}) & /@ (# + {0.1, 
            RandomVariate[NormalDistribution[0., 0.05]] + 
             0.1*f[#[[2]], r, xfunc[#[[1]] - i + 3., #[[2]], #[[3]]], 
                yfunc[#[[1]] - i + 3., #[[2]], #[[3]]]][[2]] , 
            RandomVariate[NormalDistribution[0., 0.05]] + 

             0.1*f[#[[2]], r, xfunc[#[[1]] - i + 3., #[[2]], #[[3]]], 
                yfunc[#[[1]] - i + 3., #[[2]], #[[3]]]][[
               1]] } &) , # , 10] & /@ trajectories[[-1]]]; 
   Export["~/Desktop/Rubbish/frame" <> ToString[1000 + i] <> ".gif", 
    Show[Import[
      "~/Desktop/MovieOut/frame" <> ToString[1000 + i] <> ".gif"], 
     Graphics[{Red, Disk[#, 1]} & /@ (0.9361111 # + {44., 28.} & /@ 
         trajectories[[i - 2, All, {3, 2}]]), 
      PlotRange -> {{0, 1440}, {0, 720}}]]];
   , {i, 3, 1238}], i];

Here are some frames:

enter image description here

This leads to the following animation (due to file size I only show a later part of the simulation).

enter image description here

Please have a look at the higher resolution and longer animation on Youtube.

It becomes clear that there are 5 large areas of rubbish in the oceans. This and their approximate positions are in accordance with their known positions and satellite data.

This type of simulation is very much simplified. We are making many assumptions, but by using the "observed" velocity field, we get around the fluid mechanical problems usually involved in these problems. The main idea can be used for many other problems as well. For example, we could iterate backwards to see where particles came from.

So I challenge you to simulate the following. Several fragments of the crashed Malaysian airplane have been found. Can you use the flow, invert the time direction and simulate where the parts must have come from? I wonder whether you can make better assmumptions than me (the larger fragments drift differently in the currents), and whether your point of origin correponds to mine.

I would love to hear back from you, and get ideas of how to apply this. I had a previous discussion with Vitaliy Kaurov about this, and there are certainly much nicer ways of representing the results. Any ideas are welcome.

Cheers,

Marco

POSTED BY: Marco Thiel
9 Replies

Marco and Bjoern, This is an incredible post.

About the data relevant to the MH370 flight, do you know if the data from 2014 is available? Or is there a way for you to post some small amount of data like average currents for the Indian ocean?

Searching for ocean current pictures, they seem to vary in some important details, but it appears that the wreckage appearing on Reunion island is close to several divergent integral curves. Also it looks like there is a lot of turbulence. Just based on what I see there, I'd expect sensitive dependence and the wreckage conceivably could have come from where the authorities say it crashed but just from the currents a case could be made for anywhere in the Indian Ocean, or even the South China sea which is where it was supposed to be flying.

That's just a guess, and it would be interesting to quantify that.

POSTED BY: Todd Rowland

Dear Todd,

thank you very much for your encouraging words.

About the data relevant to the MH370 flight, do you know if the data from 2014 is available? Or is there a way for you to post some small amount of data like average currents for the Indian ocean?

There is indeed some more data available on the website. If you go to this page, which I used for the download, you will see that there is data available for up to the end of October 2014.

enter image description here

You seen that the last file name ends in 20140124.nc, which indicates the date.That end data is, of course, not quite enough to follow the debris of the plane through the ocean, but you can use is for the first couple of months and then use the older data from previous years and make a sort of ensemble prediction over the flow in different years.

This is also what could be done in case of an oilspill, where flow data for predictions is not readily available. It is also interesting to study how different the paths of the particles are if you study at different times, i.e. current conditions.

Another interesting thing to try might be to use the data for different depths. If you open a file the top of it looks like this:

enter image description here

so we get all sorts of data for various different depths, not only the surface. This could be quite valuable for applications.

There are lots of other things to try. Here I have a gif-animation of the particle movements for the Fukushima scenario:

Export["~/Desktop/swarm.gif", 
 Table[Graphics[{RGBColor[i/1238., 1. - i/1238., 0.], Opacity[0.8], 
      Disk[#, 3]} & /@ trajectories[[i, 1 ;; 1000, {3, 2}]], 
   PlotRange -> {{470, 950}, {380, 570}}, ImageSize -> Full, 
   Background -> Black], {i, 1, 1236, 5}]]

enter image description here

You can use Image3D to get a representation of it, which can (after slight modification) be 3D printed.

Image3D[Import["~/Desktop/swarm.gif"]]

enter image description here

The "green top" shows the earlier times, and the "red bottom" the later ones.

Thanks a lot for your comments; I really appreciate it.

Cheers,

Marco

BTW, you can get pretty nice effects with ContourPlot:

ListContourPlot[veltottot, AspectRatio -> 1/2, ImageSize -> Full, ColorFunction -> "TemperatureMap", Frame -> None, Contours -> 7]

enter image description here

POSTED BY: Marco Thiel

Dear Marco and Bjoern,

this - as usual! - is a remarkable post! Thank you very much for sharing! I started out for this "MH370 challenge", but soon I got stuck at the rubbish problem: Assuming the plastic garbage stays next to the surface but the water moves up and down, i.e. on the surface there are (so to say) wells and sinks - or more mathematically speaking: $\nabla\vec{v}_{\scriptsize horiz}\neq0$. I thought the rubbish should accumulate on places where there is a sink. So I downloaded the data for vertical velocities for most of the year 2014. Because the data change over the year but the rubbish stay stationary, I calculated the mean value of all the data. The result looks like this (shown are only the sign of the horizontal velocity!):

enter image description here

I was very surprised to find wells and sinks mixed on a very fine scale! And the overall result is not really suited to give the nice result of the original post. Only very roughly the garbage "piles up" on places without vertical motion.

Regard -- Henrik

EDIT: There is an effect called "salt fingering", which is a mechanism for mixing water in the ocean. Maybe this is what we are seeing here. Regarding the shown horizontal motion: I expected to see the famous Humboldt current (along the west coast of South America) - what is NASA hiding?

POSTED BY: Henrik Schachner

Just wanted to say: very nice job! Perhaps you could add the border of the countries, such that there are clear distinct lines between land and sea.

POSTED BY: Sander Huisman

Dear Henrik,

Thank you very much for your post and for your time looking into this. The figure you show is really interesting. The extreme fine structure, I think I understand. It appears that the rubbish accumulate in the areas that have no pronounced fine structure (the large white areas in the oceans). I believe that there are numerical studies of much simpler situations that show that (massive?) particles often accumulate in the larger vortices, which is consistent with that I see. You are mostly seeing the very fine structure, which is the most fleeting structure. There appears to be a "fractal" pattern/scaling in the vortex size etc. It would be quite useful to see the code that you used. I

I have a rough idea what the answer to your question is, but I am not sure enough to post it here. I have a colleague, Alessandro de Moura here at the University of Aberdeen, who is an expert in this area. He has published a lot in this field, and is co-author on this paper. I have asked him for an explanation, and will post it when I get a confirmation from him.

I also fully agree with you that 3D is probably important here. It would be really nice to study the complete 3 dimensional system. As I said, the data files do contain (some of) the information we would need. I do have access to more CPU time, and one could combine efforts here in the community to sort this out.

I will try to come up with something as soon as I can.

Thanks a lot again for commenting,

Marco

POSTED BY: Marco Thiel

Dear Marco,

thank you very much for you informative answer. Of course it was naive to think there are "big" areas with wells or sinks. And one should not post a picture without any code! See the attachment, the code is quite simple. Without seeing your code I never ever would have been able to do the download - thank you for this as well! The papers you supplied seem basically to be about Hamiltonian chaos - decades ago I was quite interested in that field. I hope I find the time reading them! All in all a very nice stimulation!

Regards -- Henrik

Attachments:
POSTED BY: Henrik Schachner

Dear Henrik,

thank you very much for your notebook. It is very interesting indeed. I spoke with Alessandro, who promised to post an explanation, but apparently has not gotten around to it. I had suspected that the problem is that the flow is non-autonomous, i.e. time dependent, and that it is a 2D slice of a 3D system.

It appears that the attractors for the advection of particles can be really complicated and it is by no means clear that the "sinks" are the actual attractors. From what I understand the implications of these very complex flows are not fully understood and even much more simplified systems are rather difficult to understand.

I think that it should be possible to produce simple models to study some aspects of these flows. I am, however, not familiar with these types of simulations.

It would be great to see what type of fluid-mechanical simulations people have implemented in Mathematica. I know the demonstrations, and this beautiful 2013 blog entry.

I also had some students who used windspeed data to simulate distribution of ash from a volcanic eruption and they got excellent results. But similarly to my simulation, they basically used the flow and put particles in it.

Cheers,

Marco

POSTED BY: Marco Thiel

Dear Sander

you are quite right. I tried to do that, but had trouble aligning the images properly. It is certainly possible, but requires some fiddling with projections.

Cheers,

Marco

POSTED BY: Marco Thiel

enter image description here - this post earned "Featured Contributor" badge, congratulations !

Dear @Marco Thiel, you are the first person on Wolfram Community to be repeatedly selected as a "Featured Contributor". This is a great post and it is now featured the curated Staff Picks group. Thank you for your excellent contributions!

POSTED BY: Vitaliy Kaurov
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract