Lava? There is no way to know. Probably? If it works at all, it should work for that as well.

Let me publicly collect my thoughts. Maybe there is interest, or suggestions...

What I am going to attempt has several problems. First, the model I have implemented is for two dimensional fluids: this implies infinite depth, not infinite shallowness. Second, it is only physically accurate for microscopic scales. Third, the simulation blows up if there are high ratios between fluid densities (areas with lots of fluid, and areas with very little fluid).

Thankfully, the fluid depth should not change things too much, because I am only intending this to work under constrained conditions, and fluid depth matters less in those conditions. And we do not really care about physical accuracy, just something that looks believable that can be quickly, iteratively tweaked by an artist.

The density ratio thing is a real problem though. In this model, every cell (pixel on the heightmap) has 9 numbers associated with it, that add up to about 1. This sum represents how much total fluid there is in that cell. There is one number representing how much fluid goes into each neighboring cell (including how much fluid stays in that cell). All of these numbers are then transmitted to the neighboring cells. Each cell then updates its numbers (think of it like moving them towards a weighted average).

The problem is if one neighbor gives you a value of 10,000 and another gives you 0.0001... the update rule simply breaks. There are ways of improving this, which rely on physical theory I don't understand, but I think I can fake it. Here are some of the not-physically-accurate alterations I intend to try:

1: I will keep the simulation stable by treating the fluid density & speed differently in different contexts. The fluid update rule will see the "real" value of the fluid, but the erosion algorithm will see a (possibly nonlinear) function of that. This is to say that a cell which has a total sum of 0.8 might be treated as having no or almost no fluid in it for purposes of the erosion algorithm. A cell with 1.2 might be considered more like a lake. You can imagine that if you were to threshold the flow in this video (

https://www.youtube.com/watch?v=TkiSO3uzs48&list=PLXpzskMR1uFv51vp2FVP4KRTec9akvE8F ) that you would see river-like structures.

2: I will add a sponge layer. This slows down the computation, slightly, and uses slightly more memory, but it makes it so that fluid flows off the edge of the map with much less bounce-back. It does this by adding extra cells on the side of the simulation where the laws of physics change the further off the simulation you go. If you watch the videos I posted (

https://www.youtube.com/watch?v=3zH_zYwwdTk&index=5&list=PLXpzskMR1uFv51vp2FVP4KRTec9akvE8F ) you will notice that there is an initial shockwave, and it bounces off the side of the simulation. The shockwave was because I was lazy, and the bounceback was because I did not use a sponge layer. Neither will happen now.

3: Any regions of excessively high (or low) density fluid built up in the sponge layer will be (gently) exponentially damped. This means that excess fluid will be created or removed off the sides of the sides of the simulation as necessary to keep things stable.

4: The height of the terrain in the sponge layer is antisymmetric to the terrain perpendicular to it. That is to say, that as you move the last few pixels to the boundary you see altitudes of 108, 104, 101, 100 then the sponge layer will take on values of 99, 96, and 92.

5: I will not rely on boundaries to induce flow. I will cause flow by either (A) adding a force dependent upon the slope of the terrain at that point or (B) causing some particles to bounce back from neighboring cells that are different heights. That is, if you are at a cell/pixel of altitude 100, and you try to send some fluid to a cell of altitude 105, then some of that fluid will bounce back, and some will spill over. If you try to send very little fluid, all of it will bounce back. If you try to send fluid to a cell of altitude 95, then all of it will flow into that cell.

6: Every time fluid flows into a cell or bounces back, it erodes that cell slightly. All sediment is evenly mixed in a cell. Bounce back causes more erosion than flow. I expect to use some nonlinearity here: flow that is twice as strong is more than twice as erosive. This might be exacerbated explicitly by calculating the rotation (curl) and explicitly eroding more on one side of the map.

7: Sediment is deposited... according to some rule I have not decided on yet. Let me know if you can think of a good one. The easiest way is to have some percentage be deposited every time step.

8: I will output maps of fluid density, erosion, deposition, velocity (and maybe rotation) in addition to the heightmap itself.

9: It will take inputs of a heightmap and a precipitation map (which represents how much fluid is deposited into every cell). I will assume the precipitation map is constant valued for now: every cell gets the same rain. I could foresee more complicated inputs being possible, but I am not in a rush to implement them.

10: I do not know how I will solve the problem of mountain-tops. These areas could naturally develop fluid densities of 0, as all the fluid would flow away and little if any would return from downhill. I anticipate that the easiest thing to do will be to add a constant amount of fluid to every cell every time step, and then multiply every density by a number a little less than 1. (or decreasing it as a function of how much more than 1 it is) This should prevent any cell from getting values too big, or values too small.

Let me be very clear: I do not know if this is going to work, and I am only attempting to alter this one model because it just happens to be the one simple fluid dynamics algorithm I know. This is just an outline of my idea. There is certainly a more clever way to approach this, but I do not know what that is.

Now that I think about it, I can imagine a way to simulate a lava source... that would be interesting.