A MATLAB package for modelling geologic scale continental deformation, using the thin viscous sheet approximation

Continental collision produces some of the most striking tectonic features on the Earth's surface. For example, collision of the Indian tectonic plate into Asia at 60 Ma caused the formation of the Himalayan range and surrounding Tibetan Plateau, the highest elevation on Earth. These tectonic processes caused huge amounts of uplift and crustal thickening over a far horizontal extent. To better understand these processes, numerical models can be used to test mechanisms and physical properties.

Here, I present MATLAB package for modelling geologic scale continental deformation, using the thin viscous sheet approximation of England and McKenzie (1982). Appropriate corrections and changes are made according to later publications.

The idea behind the thin viscous sheet approximation is to represent the crust as a thin layer that is isostatically in equilibrium with the lithosphere, and that velocities within the crust do not vary with depth. This flow is incompressible and is driven by tractions and the gradients in crustal thickness. Vertical gradients of deviatoric stresses are also neglected. Considering vertical gradients in forces and rheology, it can be shown that pressure in a column of crust is

For a full description of symbols used, please refer to England and McKenzie (1982). Consider a level at the base of the lithosphere for isostatic balance (z=0). The vertical average of this pressure through the lithosphere and any topography gives

The term is calculated through the constitutive equation

This describes a power-law rheology, where n=1 is a Newtonian fluid. The parameter B is a constant which encompasses depth averaged viscosity, and is the second invariant of the strain rate tensor.

Substituting and into force-balance equations in the horizontal directions gives

In this representation, only horizontal derivatives are considered. The terms are non-dimensional and the Argand number (Ar) is defined as:

This number describes the relative contributions of forces arising from variations in crustal thicknesses and the force required to deform a fluid. For large Ar, the crustal gravitational forces succumb to deformation—the lithosphere is weak (a cheese analogue might be a ripe brie). Small Ar flows resist deformation and changes in crustal thickness, forming rigid blocks (think feta).

To constrain time-dependence, the continuity equation is written in the form:

Equations are solved on a NxN grid (e.g. 32x32), through numerical approximations to the equations for velocity and crustal thickness. The derivative terms in the RHS of the velocity equation are approximated through centre difference schemes. To solve for the velocity field on the LHS, the interior nodes of the RHS are initially set to zero and the new velocity field is inverted for using a Poisson equation solving routine.

Boundary conditions are applied to the velocity fields and crustal thickness (see above figure, from England and McKenzie (1982)). These new velocity fields are then used in the centre difference schemes to calculate a new RHS, which is solved again to attain a new solution for the velocities (again applying boundary conditions). This new solution is applied to the old solution until the solution converges. For calculations, very small convergence parameters must be used.

To advance the solution in time, as well as calculate changes in crustal thickness, numerical approximations to the continuity equation are used. The LHS is determined using forward differences, and the RHS is found using an 'upwind' scheme, as the form of the equation lends similarities to a material derivative. The timestep is chosen to satisfy the Courant-Friedrichs-Lewy criterion.

This package has been updated and benchmarked since I first wrote it in 2017, and is now ~10 times faster. The figure below shows time taken to run a 0.5 million year Newtonian simulation on a 16x16 grid, over a range of Argand numbers. Left shows the old implementation, and right shows the current. Grey lines are individual samples, and red lines are medians.

This example follows the script main.m. Further information on the functions/classes/options can be found in the headers of each file.

The constructor of class TVSfield is used to set up an object with default parameters of the simulation.

% Instantiate field object
thinViscousSheet = TVSfield();

Default settings can be changed by accessing the properties of the object. These changes must be changed before creating the grids and stencils.

thinViscousSheet.simSettings.Nx = 32;
thinViscousSheet.simSettings.Ar = 10;

Solver grids and stencils are initialized using the setupGrids() method.

% Setup grids and stencils
thinViscousSheet.setupGrids();

Specify the number of time-steps to solve for, then pass this as an argument to timeSolve(). The properties of the TVSfield object are automatically updated.

% Solve
nTimeSteps = 100;
thinViscousSheet.timeSolve(nTimeSteps);

Use the methods plot3D() and plot6() to view the results. Methods plot3D() mirrors the domain about the line x=0 for plotting purposes.

figure
thinViscousSheet.plot3D('default');

To illustrate the uses of this package I show the results of two simulations: one at Ar=1; and another at Ar=10. Both simulations are Newtonian (n=1) and are solved over a range of 5 million years. Plots of topography from plot3D() are shown below (left and right are Ar=1 and Ar=10, respectively). Note that the crustal thickening in the low Argand number case is greater but concentrated in a smaller area, in comparison to the high Argand number case.

Plots of diagnostic parameters from plot6() are shown below for the Ar=1 model. Note the changes in the strain field in the X direction and the resulting gradients crustal thickness.

Below are animations comparing the time evolution of the previous Ar=1 and Ar=10 models (left and right are Ar=1 and Ar=10, respectively). Animations are made the using the function mainAnimation() in main.m.