As promised in my last post, here is an overview of the basic code structure and interesting components that go into my planet terrain simulation.

But first, a video of my most recent run with a few examples of the successes and hiccups I will discuss below.

Basic Structure

Round: A step forward in time for which discrete physical components (tectonics, erosion, ect.) are run.

The World: Aware of all components required to progress the world forward in time. Most everything is done at this level. Moves attributes (rock, water, momentum, moisture) between plates and cells.

The Generator: Sets up the initial starting plate. Separate from the World for greater composability.

Simulation Runner: Runs the generator and world forward in time. Responsible for determining how and when to emit data between rounds.

Plates: Responsible for movement of landmasses via rotation through world space.

Plate Cells: Responsible for keeping track of columns of rock, and vertical position within the lithosphere.

World Cells: Associates plate cells with positions in world space. Determines attributes of the position based on its and its neighbors’ contents.

The Grid: Contains vertex data as well as information on which vertices are neighbors. Each cell is associated with a particular vertex and rotated space.

Time Complexity

The time complexity of the current build is O(n*sqrt(n)). Each part of the simulation is at most O(n), and the time step used for each round scales with the radius of simulation cells as sqrt(n).

Grid Creation

The grid is currently constructed with this library. Starting with an icosahedron each face is divided into a large number of equilateral triangles and projected onto a sphere. This projection results in a slightly larger edge length for triangles originating near the center of a face of the icosahedron. The resulting artifacts appear to be minimal, but should be addressed eventually.

I originally created this library for creating spherical hex based game grids.

Meteor Initialization

The current generator initializes the world in the style of the late heavy bombardment period on Earth. Thousands of craters are created at random points with radii ranging from 200–2500km. As of my first implementation, mass is ejected within the radius such that depth equal to 1/20 of the radius is removed near the center, decreasing in a quadratic manner to zero at the outer rim. The total ejected mass is then deposited according to a Gaussian distribution centered along the outer rim (this should be replaced with a gamma distribution or the like, so that the crater rim is not a giant cliff). To maintain a constant crust mass, the weight of each receiving cell is first calculated and then normalized to the total weight before ejected mass is deposited.

Craters are currently implemented as additive elevation changes, so newer craters can be seen overlapping instead of occluding old ones. Eventually later impacts should take into account existing craters. Also a better crater shape needs to be implemented, but this model works for now.

Plate Movement

Plate movement is handled via a rotation matrix, angular speed, and rotation pole. The speed and pole could be combined, but I prefer the clarity of separate values as no optimization is needed here.

Plate Casting

At the beginning of each round, every plate cell is ‘cast’ to the nearest world cell. The new nearest cell is calculated from the set of the previous nearest cell and its neighbors. If this results in a different nearest cell the process is repeated. Distance squared is used for comparison to avoid unnecessary square roots.

With some preprocessing of the vertex grid, I should be able to determine a distance within which a cell is know to be closest. This distance would allow the above algorithm to terminate early if a world cell was found within that distance.

Plate Boundaries

A giant mess. World cells without any corresponding rock from a plate cell are filled if possible, and world cells with rock from multiple plates transfer rock from one plate to another based on plate index. I will write more about this process once I have fixed it.

World Level Smoothing

Occurs before segments of a round that require full coverage of world cell attributes (e.g. elevation for erosion). For elevation it is a simple recursive algorithm which replaces invalid world cell elevations with the average elevation of the cell’s valid neighbors.

Rock Erosion

Erosion happens in two parts. In the first, each world cell creates sediment in neighboring cells based linearly on height difference. Precipitation dependence will be added once I have a model for it.

The second is a temporary measure to control mountain height. As a cubic function of height, each cell erodes its rock into sediment such that at 4,000 meters 1/400 of the elevation is eroded per million years, and up to a cap of 1/2 of a cell’s elevation at ~26,000 meters.

Sediment Transport

Once sediment has been created, an acyclic directed graph is constructed such that each cell above the continental shelf (currently static at 300 meters below sea level) passes incoming sediment and a portion of its own sediment to downhill cells.

The amount of sediment taken from each cell is based on the steepest slope with downhill neighbors such that sediment moved matches the volume of suspended sediment found in corresponding areas along the Mississippi River basin. On first pass, the closest equation that fits this data and behaves well at higher gradients is as follows:

Remove Sediment(meters) = 15237.4*ln(slope+1) + 170378.0*ln(slope+1)²

Other rivers on Earth will be added to the model once precipitation is more than a world constant.

The fraction of sediment moved to each downhill node is weighted by the squared difference in elevation. The idea is that the stability of river paths depends on the variance in slope to downhill nodes. While I have no research to back this model, it produces nicer looking results than the linear model I have used previously.

As I haven’t added sediment deposition yet, everything ends up in large piles at the end nodes of the graph. As a result, ‘marching’ mountains can be seen as these piles move towards higher elevations.

How the mountains march.

Homeostasis

Here I determine the actual elevation of each plate cell by floating each cell in the slightly denser mantle. Cells do not affect their neighbors, which might be the cause of some of my elevation problems. I’m looking at modifying this phase towards an elastic lithosphere model.

Angular Momentum

A conservative angular momentum model is in place for each plate. When rock moves from one plate to another, so does its momentum. Additional collision momentum needs to be transferred such that continents eventually dock, however I have not yet found a good method to determine realistic rates for this transfer.

During phases that transfer momentum, the total transfer from a source plate to a destination plate is summed as a mass*radius scalar and multiplied by the relevant angular momentum vector at very the end.

What’s Next?

A new tectonic model is in the works. After I complete that, sediment deposition will be added.

A brief overview of the scientific papers I have used so far is still coming. If you find this project interesting, please follow! I will also be posting short articles as I update and add to the existing models and algorithms.