Elevations blow up!

[astrographer]

Acholas-599Ma.zip As shown in the attached images, the elevations are getting to excessive levels. In this example, the minimum is -74.4km and the maximum is 21.7km, clearly unrealistic and it can spoil any attempt at a ~realistic~… er… plausible climate model.

I'm not sure why it is doing this, so I'm not able to offer any means of fixing it backed by rigorous research. That said, although it lacks any backing, perhaps the simulator could periodically rescale elevations back to a -11000 to 8500m range, and downscaling crust thicknesses commensurately. The problem seems to build up slowly, so perhaps every 10-100 cycles.

Clearly, a much better solution would be to track down why the values are getting so extreme…

[astrographer]

In terms of value, I'd put this issue at the top of my list, but it feels like it may be a hard fix…

[davidson16807]

That is most disconcerting.

OK, so elevation is simply displacement minus sealevel, where sealevel remains constant. Displacement is a function of the thickness and density of the crust reflecting isostasy, so dense thin crust ought to be at lowest elevation. Oceanic and continental crust behave differently in terms of density and thickness. Oceanic crust density depends on age, so older crust is more dense. Continental crust never changes density, so it's age is never really tracked. Crust thickness only changes when plates overlap. When old ocean crust overlaps, it is destroyed. Continental crust is a conserved quantity, so it should never be created or destroyed.

It could be some sort of subtlety concerning age. "Crust age" should always refer to the age of oceanic crust, but there might have been a slip up at some point. It's also possible the reset that occurs at the start of a new supercontinent cycle is changing age where new crust overlaps old crust, so old crust becomes new or new crust becomes old. These are just guesses. I'll have to look into this.

[davidson16807]

OK, I was able to reproduce with the seed

tectonics.js/index.html?seed=Acholas

but this behavior is a lot more sporadic than I expected. It doesn't always blow up at 600My, and even when it does, it still has the chance of returning to a normal state.

Also, density is really wierd when this happens, higher than 9000 kg/m^3. I didn't think it was possible for the model to achieve that.

I wonder if this is more something to do with the display, or maybe it's some nonsense that results from sharing rasters for performance optimization. I'm starting to think its not something to do with the model itself.

I wonder if you or I could generate a save file of when the problem occurs.

[davidson16807]

Another clue: I reproduced the problem again, and this time I tried generating a save file. However, when I tried to load up that save file the display went back to normal. So it's not a problem with model state. I think this is a problem with the display.

[davidson16807]

Nope, nevermind. I setup a debug flag to trigger when elevation gets lower than -10km. The offending grid cell I found has a thickness of -51km: of which 7km is sima, and -58km is sial. That's bananas.

So the problem is sial can reach negative numbers. That's doubly strange because we never subtract anything from it, like, anywhere.

[astrographer]

Just to be clear, sima is ocean floor/deep crust basalt and sial is continental crust, correct?

On Dec 15, 2017, at 10:55 AM, Carl Davidson notifications@github.com wrote:

Nope, nevermind. I setup a debug flag to trigger when elevation gets lower than -10km. The offending grid cell I found has a thickness of -51km: of which 7km is sima, and -58km is sial. That's bananas.

So the problem is sial can reach negative numbers, which is doubly strange because we never subtract anything from it anywhere.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/davidson16807/tectonics.js/issues/17#issuecomment-352083938, or mute the thread https://github.com/notifications/unsubscribe-auth/AGayL1WIqT9DxIk-9N5Pn1DYI41VgKOhks5tAsChgaJpZM4RCgRG.

[davidson16807]

Correct, continental. Actually, I mispoke - we do subtract it in places (I was confusing it with ocean crust)

I set up a custom display to highlight the negative regions:

view.setScalarDisplay( new ScalarHeatDisplay( { min: '-1.', max: '0.',  
        getField: function (crust) {
            return crust.sial;
        } 
    } ) 
);

[davidson16807]

Oh, ok. This must be an issue with erosion. Erosion subtracts sial from high elevations and adds it to low elevations. But if a high elevation exists with no sial, then the sial gets pushed into the negatives.

This is consistent with the fact that all regions with negative sial are in the ocean - the problem only occurs when ocean crust is higher than neighboring crust.

This also explains why the issue corrects itself over time - an area with negative sial becomes very low in elevation, so it fills back up with sial at first opportunity. And despite the wierdness, sial remains a conserved quantity, so it doesn't invalidate the model in that regard.

There needs to be a more robust way to track a conserved quantity like sial throughout the model, sort of like accounting software.

[davidson16807]

Status update:

The problem goes away when I disable erosion. This confirms my suspicion that it's something to do with erosion.

I started a new branch, called "sediment". I re-read Simoes 2010 and came to conclude there's a lot wrong with the current implementation of erosion. I'm not even sure I can enumerate it all:

• erosion occurs on ocean floor, despite there already being code to prevent that from happening

• The model tracks rock by thickness, not mass. This makes it much harder to guarantee mass conservation, and makes it harder to expand out to include other rock types, since other rock types have their own unique densities.

• Only one part of the Simoes model was implemented, and this part actually doesn't express mass transport at all. It actually expresses how rock weathers down to sediment. There's a separate component that governs how sediment is transported, and this wasn't implemented. My implementation just assumes all sediment is immediately transported to adjacent cells, where it becomes rock again. This meant the entire model could simplify down to diffusion. Come to think of it, I'm surprised this model worked as well as it did for so long. Which leads me to the next problem...

• Mass transport is modeled by the diffusion equation. This is fine as long as there's only one rock type, but it breaks down once you add multiple rock types like sima/sial. The problem we observed is a direct consequence of this.

Rather than solve these individually I'm ripping out the erosion model and re-implementing it. It's going to be much closer to the way Simoes described it, except for two differences:

• it's going to track mass, not thickness

• it's going to be expressed in a way that's more in line with the diffusion-convection equation.

Tracking mass means there's less that can go wrong during conversion from bedrock to sediment. Expressing the model in terms of the diffusion-convection equation means it's going to resemble a lot of other models that I foresee needing in the future, so there's more opportunity for code reuse.

A side-effect of this means I'm going to be implementing sediment from issue #8.

[davidson16807]

Another status update:

The updated erosion model is implemented on the sediment branch, but the results are even more unstable than the original model. There's a lot of underlying functionality that's never been tested before. I don't think I've ever properly tested the functions for calculating gradient and divergence. I've already discovered several bugs in those functions, and I suspect there are more. I actually think I'll need to completely rewrite the gradient calculation. If anyone wants to offer a performant implementation for calculating gradients on an unstructured mesh, please do.

[davidson16807]

I did some research on calculating gradients from unstructured meshes. There's a pretty popular implementation that uses the Gauss-Green theorem, and I'm trying it out. Results so far are promising. Sial/sima retains positive values throughout simulation. Sediment is still reaching into the negatives, but the model is stable, and it never seems to cause thickness or density to reach into the negatives.

[davidson16807]

The implementation is very stable when plate motion is 0 but blows up whenever I try to turn plate motion back on. Also, erosion appears much slower, now, almost ineffective. Hrm.

[davidson16807]

I've taken the pragmatic approach. I'm starting from scratch on a new branch. There were too many sweeping changes on the last branch, I couldn't make sense of what's broken anymore. This time, I'm addressing the immediate problem using the existing approach, but slightly modified to prevent negative values. If a cell's outbound sediment transfer exceeds the quantity of sial available, the outbound sediment transfer scales back to match the supply.

Now, this should absolutely work in theory. However, I'm still getting negative values soon after startup. Here's the kicker: the negative values go away immediately after turning plate motion off. So the problem goes away if erosion is turned off, and the problem goes away if plate motion is turned off. This is a problem having to do with how the submodels work together. I'm not sure yet how this comes to be the case.

[davidson16807]

OK, it's nothing to do with how deltas are calculated. Erosion is calculated once plate motion is already applied, like it should. I think the remaining negative values have something to do with the way erosion is applied to each plate.

Erosion is calculated once for the entire globe, because it really is a global thing - sediment ought to be transported across plate boundaries. However, sial/sima is represented individually for each plate. So for each cell the model applies erosion to whatever is the topmost plate. But the function used to lookup the topmost plate is imprecise. Sometimes near plate boundaries it gets the wrong plate. In that case, it could take away more sediment than what's already there. This is my best guess.

[davidson16807]

OK, I think I understand the problem well enough, and I believe the remaining effect only applies to cells that are not associated with any plate, i.e. it's pretty negligible as far as what the model says about the actual planet. I'm fine with clamping erosion after its calculation so that it doesn't exceed available sial. Fix is out on the sediment2 branch, and it appears sial never reaches into the negatives. Results appear very stable.

Give it a try if you'd like. This should be enough to close the issue, but it's put me in a head space where I can address some other enhancements related to erosion. I still want to:

• represent sial/sima in the model using mass, not thickness

• represent sediment within the model as per issue #8

• represent weathering and erosion separately, as per the Simoes model that I originally intended to use, or possibly switch to using the convection-diffusion equation. If I can get the convection-diffusion equation to work it opens up new vistas for modeling the concentration of specific types of minerals and maybe even the spread of biological taxons. However, this might get put on the backburner until I can improve climate modeling and biome distribution.

[davidson16807]

Aside: this isn't necessary to close out the story, but I did some work to reparameterize the model. Previously, some elevations reported by the chart were just absurd, like, two or three times the height of Everest. The model looks much more Earthlike, now. Changes are out on master.

[astrographer]

Sweet! I'll check on it ASAP.

On Feb 9, 2018, at 3:33 PM, Carl Davidson notifications@github.com wrote:

Aside: this isn't necessary to close out the story, but I did some work to reparameterize the model. Previously, some elevations reported by the chart were just absurd, like, two or three times the height of Everest. Results look much more Earthlike, in general. Changes are out on master right now.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or mute the thread.

[davidson16807]

Hrm, I don't think sial is being conserved now. I think the model is popping it into existence somewhere. Elevations tend to reach some very unrealistic values if you run the simulation for long enough, again way past Everest. The conservation violation only occurs when both erosion and plate motion is turned on.

[davidson16807]

I pushed a commit to ensure conservation. Elevations never get much past 20km, or twice the height of Everest. That's still implausibly high, but it's nothing the model didn't do previously. Nuts to this, I'm working on sediment.

[davidson16807]

I'll close this out if there are no objections.