Notebook - Mandyoc tutorial


Mandyoc: MANtle DYnamics simulatOr Code

Mandyoc is a finite element code written on top of the PETSc library to simulate thermo-chemical convection of the Earth’s mantle. Different linear and non-linear rheologies can be adopted, appropriately simulating the strain and stress pattern in the Earth’s crust and mantle, both in extensional or collisional tectonics. Additionally, the code allows variations of boundary condition for the velocity field in space and time, simulating different pulses of tectonism in the same numerical scenario

To simulate mantle thermochemical convection, we adopted the formulation for non-Newtonian fluids together with the Boussinesq approximation to solve the following equations of conservation of mass, momentum and energy:

ui,i=0u_{i,i} = 0
σij,j+giρ0(1α(TT0))=0\sigma_{i j,j} + g_{i} \rho_{0}(1 - \alpha(T - T_{0})) = 0
Tt+uiT,i=κT,ii+HcpαTguecp\frac{\partial T}{\partial t} + u_{i} T_{,i} = \kappa T_{,ii} + \frac{H}{c_p} - \frac{\alpha T g u_{e}}{c_{p}}

where uu is the velocity in the ii direction, gg is the gravity acceleration, ρ0\rho_{0} is the reference rock density at temperature T0T_0, α\alpha is the coefficient of volumetric expansion, TT is the temperature, κ\kappa is the thermal diffusivity, HH is the rate of radiogenic heat per unit of mass, cpc_{p} is the specific heat, δij\delta_{ij} is the Kronecker delta, and σij\sigma_{ij} is the stress tensor:

σij=Pδij+η(ui,j+uj,i)\sigma_{i j} = -P \delta_{i j} + \eta (u_{i,j} + u_{j,i})

where PP is the dynamic pressure and η\eta is the effective viscosity of the rock.

The equations of conservation of mass, momentum and energy are solved using the finite element method. Mandyoc uses quadrilateral elements for two dimensional grids. The mass and momentum equations subsection presents the numerical methods used to solve the mass and momentum equations and the energy equation subsection shows how an implicit formulation was used to solve the energy equation.

To simulate any scenario, the user must provide the parameter file param.txt and, if set, the necessary ASCII files with the initial temperature field, velocity field and/or the initial interfaces of the model. The flowchart showing the steps Mandyoc takes to solve the equations of conservation of mass, momentum and energy.

It shows that once the code starts running and the input files are read, Mandyoc uses the effective viscosity field η\eta to calculate the velocity field uu and checks if the convergence condition satisfies the tolerance toltol.

While the minimum tolerance is not reached, Mandyoc utilizes the Uzawa's method to iteratively calculate a new uu and PP fields. The updated fields modify the viscosity field η\eta, which in turn disturbs the velocity field again. These fields are updated until tolerance is reached.

Additionally, the compositional factor CC is evaluated for an advection as:

Ct+uiC,i=0\frac{\partial C}{\partial t} + u_{i}C_{,i} = 0

Its solution is calculated placing randomly a number of particles within each finite element of the mesh, which are displaced based on the adjacent node velocity values. The individual value for each particle is obtained by linear interpolation of the node values.

List of scientific works done with Mandyoc or its earlier version:

Showcase: Continental rift

This example simulates the evolution of divergent margins, taking into account the plastic rheology and the sin-rift geodynamics. The idea of this tutorial is to explain how to generate the initial model to reach the following result:

Input files

There are a few ways to run this model. The one that will present comes from the assumption that the geometry, temperature and velocity at the beginning of the simulation are know in a reasonable fashion. Therefore, we must provide them as input files, which are:

  • an interface file with the boundary of the different layers and their composition information,
  • an initial temperature file with the initial temperature field,
  • and an initial velocity file with the initial velocity field

These files describe the following setup model for the continental rift:

# Python packages:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

# Import wornings package to ignore any wrnings:
import warnings

# Import some functions that I created for this tutorial:
from functions import (

The domain of the model comprises 1600 x 300 km2 that is discretized by elements of 1 x 1 km2:

lx, lz = 1600.0e3, 300.0e3

nx, nz = 1601, 301

The regular meshgrid can be easily created using:

x = np.linspace(0, lx, nx)
z = np.linspace(-lz, 0, nz)

xx, zz = np.meshgrid(x, z)

Interface file

For this simulation, the domain is subdivided into a few lithologial units:

  • asthenosphere,
  • lithospheric mantle,
  • weak seed,
  • lower crust,
  • upper crust, and
  • sticky air.

Each of them with their different compositional properties.

Additionally, to ensure the nucleation of rifting at the center of the numerical domain, a weak seed will add in the lithospheric mantle.

h_sair = 40.0e3

h_lower_crust = 20.0e3

h_upper_crust = 20.0e3

h_litho = 130.0e3

# Seed depth position under the crust: 
z_seed = 13.0e3

The different lithologies can be expressed in a dictionary:

interfaces = {
    "litho": -np.ones(nx) * h_litho,  
    "seed_base": -np.ones(nx) * (z_seed + h_lower_crust + h_upper_crust),
    "seed_top": -np.ones(nx) * (z_seed + h_lower_crust + h_upper_crust),
    "lower_crust": -np.ones(nx) * (h_lower_crust + h_upper_crust),
    "upper_crust": -np.ones(nx) * h_upper_crust,
    "air": np.zeros(nx),      
for i in interfaces:
    interfaces[i] = interfaces[i] - h_sair

The considerations above result in the geometry shown in the following plot:

fig, ax = plt.subplots(figsize=(16, 8))

for label, layer in interfaces.items():
        ax.plot(x/1e3, (layer + h_sair)/1e3, label=f"{label}")

ax.set_yticks(np.arange((-lz + h_sair)/1e3, h_sair/1e3, 10))
ax.set_xlim([0, lx/1e3])
ax.set_ylim([(-lz + h_sair)/1e3, h_sair/1e3])
ax.set_xlabel("x [km]")
ax.set_ylabel("Depth [km]")

plt.legend(loc='upper right')
<Figure size 1152x576 with 1 Axes>

The weak seed can be created:

# Seed horizontal position:
x_seed = 750.0e3
# Seed thickness:
h_seed = 6.0e3
d_x = lx / (nx - 1)

cond = (x_seed - d_x <= x) & (x < x_seed + d_x)

# Add a thickness in the seed layer:
interfaces["seed_base"][cond] = interfaces["seed_base"][cond] - h_seed//2

interfaces["seed_top"][cond] = interfaces["seed_top"][cond] + h_seed//2

Plot again the interfaces:

fig, ax = plt.subplots(figsize=(16, 8))

for label, layer in interfaces.items():
    ax.plot(x/1e3, (layer + h_sair)/1e3, label=f"{label}")

ax.set_yticks(np.arange((-lz + h_sair)/1e3, h_sair/1e3, 10))
ax.set_xlim([0, lx/1e3])
ax.set_ylim([(-lz + h_sair)/1e3, h_sair/1e3])
ax.set_xlabel("x [km]")
ax.set_ylabel("Depth [km]")

plt.legend(loc='upper right')
<Figure size 1152x576 with 1 Axes>

Now we can create the interface file. This file contains the layers properties and the interface's depth between these layers.

We are going to adopt a visco-plastic rheology where the effective viscosity η\eta combines nonlinear power law viscous rheology and a plastic yield criterion. The viscous component is given by:

ηvisc=CA1/nε˙(1n)/nexpQ+VPnRT\eta_{visc} = C A^{-1/n} \dot{\varepsilon}^{(1-n)/n} \exp{\frac{Q + V P}{n R T}}

So, the layers properties that we need to specify are:

  • Compositional factor (CC)
  • Density (rhorho)
  • Radiogenic heat (HH)
  • Pre-exponential scale factor (AA)
  • Power law exponent (nn)
  • Activation energy (QQ)
  • Activation volume (VV)

To simulate the free surface, the sticky air approach is adopted, taking into account a 40 km thick layer with a relatively low viscosity material but with a compatible density with the atmospheric air.

AsthenosphereLitho mantleSeedLitho mantle over the seedLower crustUpper crustAir
# Define the density and the radiogenic heat for the upper and lower crust:
rho_uc = 2700.0
rho_lc = 2800.0
Huc = 2.5e-6 / rho_uc
Hlc = 0.8e-6 / rho_lc
# Define the layer properties:
layer_properties = f""" 
        C   1.0       1.0        0.1        1.0        1.0         1.0         1.0
        rho 3378.0    3354.0     3354.0     3354.0     {rho_lc}    {rho_uc}    1.0
        H   0.0       9.0e-12    9.0e-12    9.0e-12    {Hlc}       {Huc}       0.0
        A   1.393e-14 2.4168e-15 2.4168e-15 2.4168e-15 8.574e-28   8.574e-28   1.0e-18
        n   3.0       3.5        3.5        3.5        4.0         4.0         1.0
        Q   429.0e3   540.0e3    540.0e3    540.0e3    222.0e3     222.0e3     0.0
        V   15.0e-6   25.0e-6    25.0e-6    25.0e-6    0.0         0.0         0.0
# Create and save the interface file:
with open("interfaces.txt", "w") as f:
    # Write the layer properties
    for line in layer_properties.split("\n"):
        line = line.strip()
        if len(line):
            f.write(" ".join(line.split()) + "\n")

    # Write the interfaces of the layers adding the thickness of the sticky air layer:
    data = np.array(tuple(interfaces.values())).T
    np.savetxt(f, data, fmt="%.1f")

Initial temperature file

The initial temperature structure is depth dependent and is 0°C at the surface and 1300°C at the base of the lithosphere at 130 km. With these boundary conditions, the initial temperature structure in the interior of the lithosphere is given by the solution of 1-D diffusion equation:

κ2T(z)z2+H(z)cp=0\kappa \frac{\partial^2 T(z)}{\partial z^2} + \frac{H(z)}{c_p} = 0

where H(z)H(z) is the internal heat production of the different layers.

The sublithospheric temperature follows an adiabatic increase up to the bottom of the model:

T=TpexpgαzcpT = T_p \exp{\frac{g \alpha z}{c_p}}

Where TpT_p is the potential temperature for the mantle, gg is the gravity value, α\alpha is the volumetric expansion coefficient, cpc_p is the specific heat capacity.

  • κ=1x106\kappa = 1x10^{-6} m2/sm²/s
  • cp=1250c_p = 1250 Jkg1K1J kg^{-1} K^{-1}
  • Tp=1262T_p = 1262 °C
  • g=10g = -10 m/s2m/s²
  • α=3.28x105\alpha = 3.28x10^{-5} K1K^{-1}
kappa = 1.0e-6
ccapacity = 1250
tem_p = 1262
g = -10
alpha = 3.28e-5
# Temperature when z is over 130 km:
temp_z = 1300 * (-z - h_sair) / (h_litho)

# Sublithospheric temperature:
temp_adiabatic = tem_p / np.exp(g * alpha * (-z - h_sair) / ccapacity)

# T = 0 when z is over zero:
temp_z[z + h_sair >= 0.0] = 0.0

# When z is below 130 km (Lithosperic depth), the T = temp_adiabatic:
temp_z[temp_z > temp_adiabatic] = temp_adiabatic[temp_z > temp_adiabatic]

Now, we will apply the 1-D diffusion equation in the model.

First, create the internal heat production model:

# Initialize H as a grid of zeros:
H = np.zeros_like(temp_z)

# Add the H value for the upper crust:
cond = (z <= interfaces["air"][0]) & (z > interfaces["upper_crust"][0])
H[cond] = Huc

# Add the H value for the lower crust:
cond = (z <= interfaces["upper_crust"][0]) & (z > interfaces["lower_crust"][0])
H[cond] = Hlc

Apply the diffusion equation for 10000 years:

dt = 10000
dz = lz / (nz - 1)

# Apply the 1-D diffusion equation:
temp_z = diffusion_equation(

Save the temperature field in a file:

temp_z = np.ones_like(xx) * temp_z[:, None]


# Save the initial temperature file
    np.reshape(temp_z, (nx * nz)), 
(301, 1601)

Plot the temperature field and a thermal profile:

fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(16, 8), sharey=True)

# Plot temperature field:
img_temp = axs[0].contourf(
    (zz + h_sair)/1.0e3,
    levels=np.arange(0, 1610, 100)
cbar = fig.colorbar(img_temp, orientation='vertical', ax=axs[0])
cbar.set_label("Temperature [°C]")

axs[0].set_ylim([(-lz + h_sair)/1e3, h_sair/1e3])
axs[0].set_ylabel("Depth [km]")
axs[0].set_xlabel("x [km]")

# Plot temperature profile:
axs[1].set_title("Thermal profile [°C]")
axs[1].plot(temp_z[:, 0], (z + h_sair)/1.0e3, "-k")
# Add interfaces:
for label, layer in interfaces.items():
    axs[1].plot(x/1e3, (layer + h_sair)/1e3, label=f"{label}")

axs[1].set_ylim([(-lz + h_sair)/1e3, h_sair/1e3])
axs[1].set_xlabel("Temperature [°C]")
axs[1].set_ylabel("Depth [km]")
plt.legend(loc="lower left")
<Figure size 1152x576 with 3 Axes>

Initial velocity field

The velocity field simulate the lithospheric stretching assuming a reference frame fixed on the lithospheric plate on the left side of the model, and the plate on the right side moves rightward with a velocity of vrift=1v_{rift} = 1 cm/year. The velocity field in the left and right boundaries of the model is chosen to ensure conservation of mass and is symmetrical if the adopted reference frame movies to the right with a velocity of vrift/2v_{rift}/2 cm/year relative to the left plate. The inflow on the vertical boundaries below the lithosphere compensates for outflow.

So, the horizontal velocity field along the left and right borders of the domain are:

velocityx(x=0)={0,z>hchsairvz+hc+hsairha,zhchsairvelocity_x(x = 0) = \left\{ \begin{array}{r l} 0, & \text{$z > -h_c - h_{sair}$} \\ -v \frac{z + h_c + h_{sair}}{h_a}, & \text{$z \leq -h_c - h_{sair}$} \end{array} \right.
velocityx(x=lx)={0,z>hsairvrift,z>hchsairvz+hc+hsairha+vrift,zhchsairvelocity_x(x = lx) = \left\{ \begin{array}{rl} 0, & \text{$z > h_{sair}$} \\ v_{rift}, & \text{$z > -h_c - h_{sair}$} \\ v \frac{z + h_c + h_{sair}}{h_a} + v_{rift}, & \text{$z \leq -h_c - h_{sair}$} \end{array} \right.

where v=vriftlzhsairhav = v_{rift} \frac{lz - h_{sair}}{h_a} so that:

lz0velocityx(x=0)=lz0velocityx(x=lx)=vha2 \int^{0}_{-lz} velocity_x(x = 0) = \int^{0}_{-lz} velocity_x(x = lx) = \frac{v h_a}{2}

hc=150h_c = 150 km is the thickness of the upper layer with constant velocity, corresponding to the lithosphere hlitho=130h_{litho} = 130 km and 20 km of the asthenosphere and ha=lzhchsairh_a = lz - h_c - h_{sair}

# Converting 1 cm/year to m/s:
v_rift = 1e-2 / (365 * 24 * 3600)
h_c= 150.0e3
h_a = lz - h_sair - h_c
v = v_rift * (lz - h_sair) / h_a
# Creating horizontal and vertical velocity:
velocity_xx, velocity_zz = np.zeros_like(xx), np.zeros_like(xx)

Velocity for the left side (x == 0):

cond = (xx == 0) & (zz <= -h_c - h_sair)
velocity_xx[cond] += -v * (zz[cond] + h_c + h_sair) / h_a

Velocity for the rigth side (x == lx):

cond = (xx == lx) & (zz > -h_c - h_sair)
velocity_xx[cond] += v_rift

cond = (xx == lx)& (zz <= -h_c - h_sair)
velocity_xx[cond] += v * (zz[cond] + h_c + h_sair) / h_a + v_rift

# Consider 10 km of the air layer with constant velocity:
frac_air = 10e3
velocity_xx[zz >= -h_sair + frac_air] = 0

Due to the mass conservation is assumed, the sum of the integrals over the boundaries (material flow) must be zero.

# For the left side:
v0 = velocity_xx[xx == 0]
sum_velocity_left = np.sum(v0[1:-1]) + (v0[0] + v0[-1]) / 2.0

# For the right side:
vf = velocity_xx[(xx == lx)]
sum_velocity_right = np.sum(vf[1:-1]) + (vf[0] + vf[-1]) / 2.0

diff = (sum_velocity_right - sum_velocity_left) * dz

print("Sum of the integrals over the boundary is:", diff)
Sum of the integrals over the boundary is: 3.012430238457634e-06

When the sum of the integrals on the boundaries is not zero due to rounding errors, we add the difference back on the top to compensate make it zero. Despite not being negligible, this is reasonable because it is a very small correction.

velocity_zz[zz == 0] = -diff / lx

Plot the velocity:

fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(9, 9), sharey=True)

# Plot horizontal component:
axs[0].plot(velocity_xx[:, 0], (z + h_sair)/1e3, "k-", label="left side")
axs[0].plot(velocity_xx[:, -1], (z + h_sair)/1e3, "r-", label="right side")

axs[0].set_ylim([(-lz + h_sair)/1e3 , h_sair/1e3])
axs[0].set_xlim([-8e-10, 8e-10])
axs[0].set_xlabel(" Velocity [m/s]")
axs[0].set_ylabel("Depth [km]")
axs[0].set_title("Horizontal component of velocity")

# Plot vertical component
axs[1].plot(velocity_zz[:, 0], (z + h_sair)/1e3, "k-", label="left side")
axs[1].plot(velocity_zz[:, -1], (z + h_sair)/1e3, "r-", label="right side")

axs[1].set_xlim([-8e-10, 8e-10])
axs[1].set_xlabel(" Velocity [m/s]")
axs[1].set_title("Vertical component of velocity")
<Figure size 648x648 with 2 Axes>

Create and save the initial velocity file:

velocity = np.zeros((2, nx * nz))

velocity[0, :] = np.copy(np.reshape(velocity_xx, nx * nz))
velocity[1, :] = np.copy(np.reshape(velocity_zz, nx * nz))

velocity = np.reshape(velocity.T, (np.size(velocity)))

np.savetxt("input_velocity_0.txt", velocity, header="v1\nv2\nv3\nv4")

Create the parameter file

The parameter file contains the information that is necessary for the simulation to run. There are many parameters to define, so we will only show/discuss some of the important ones for this simulation. But I invite you to visit the Mandyoc documentation website to know more about the parameter file.

  • Define the number of Lagrangian particles in each element, particles_per_element = 40.

  • Specify the maximum time and step of the simulation with step_max = 5000 and time_max = 100.0e6 (100 My).

  • For the temperature boundary conditions, the temperature values is fixed on all sides during the simulation, with the nodes of the model in the air maintained at 0°C.

  • For the velocity boundary conditions, free slip is applied on the top and bottom of the numerical domain.

  • To avoid a symmetric, laterally homogeneous model, we introduced a random perturbation of the initial strain in each finite element of the model.

  • Set the pre-defined rheology model to use during simulation. In this case: rheology_model = 9

params = f"""
# Geometry
nx                                  = {nx}          # n. of nodes in the horizontal direction
nz                                  = {nz}          # n. of nodes in the vertical direction
lx                                  = {lx}          # extent in the horizontal direction
lz                                  = {lz}          # extent in the vertical direction

# Simulation options
denok                               = 1.0e-15       # tolerance criterion for the Uzawa's scheme, default is 1.0e-4
Xi_min                              = 1.0e-7        # tolerance criterion for the convergence of the non-linear flow, default is 1.0e-14
random_initial_strain               = 0.3           # non-dimensional value for the initial strain perturbation for the entire domain, default is 0.0 
initial_dynamic_range               = True          # method to smoothen convergence of the velocity field in scenarios with wide viscosity range, default is False [True/False], Gerya (2019)
# Particles options
particles_per_element               = 40            # n. of Lagrangian particles in each element, default is 81
particles_perturb_factor            = 0.7           # indicates the amount of perturbation of the initial location of the particles relative to a regular grid distribution. Default is 0.5 [values are between 0 and 1]

# Time constrains
step_max                            = 5000          # maximum time-step of the simulation [steps]
time_max                            = 100.0e6       # maximum time of the simulation [years]
dt_max                              = 10.0e6        # maximum time between steps of the simulation [years]
step_print                          = 10            # make output files every <step_print>
sub_division_time_step              = 0.5           # re-scale value for the calculated time-step, default is 1.0

# Viscosity
viscosity_reference                 = 1.0e26        # reference mantle viscosity [Pa.s]
viscosity_max                       = 1.0e25        # maximum viscosity during simulation [Pa.s]
viscosity_min                       = 1.0e18        # minimum viscosity during simulation [Pa.s]
viscosity_per_element               = constant      # sets if viscosity is constant or linearly variable for every element, default is variable [constant/variable]
viscosity_mean_method               = arithmetic    # defines method do calculate the viscosity for each element, default is harmonic [harmonic/arithmetic]
viscosity_dependence                = pressure      # defines if viscosity depends on pressure or depth, default is depth [pressure/depth]
# External ASCII inputs/outputs
interfaces_from_ascii               = True          # set if interfaces between lithologies are read from an ASCII file (interfaces.txt), default is False [True/False]
n_interfaces                        = {len(interfaces.keys())}            # set the number of interfaces to be read from the interfaces ASCII file (interfaces.txt) 
temperature_from_ascii              = True          # set if initial temperature is read from an ASCII file (input_temperature_0.txt), default is False [True/False]
velocity_from_ascii                 = True          # set if initial velocity field is read from an ASCII file (input_velocity_0.txt), default is False [True/False]
variable_bcv                        = False         # allows velocity field re-scaling through time according to an ASCII file (scale_bcv.txt), default is False [True/False]
multi_velocity                      = False         # set if boundary velocities can change with time from ASCII file(s) (multi_veloc.txt and additional input_velocity_[X].txt files), default is False [True/False]
binary_output                       = False         # set if output is in binary format, default is False [True/False]
print_step_files                    = True          # set if the particles position are printed to an output file, default is True [True/False]

sticky_blanket_air                  = True          # default is False [True/False]
precipitation_profile_from_ascii    = False         # default is False [True/False]
climate_change_from_ascii           = False         # default is False [True/False]

# Physical parameters
temperature_difference              = 1500.         # temperature difference between the top and bottom of the model (relevant if <temperature_from_ascii> is False) [K]
thermal_expansion_coefficient       = 3.28e-5       # value for the coefficient of thermal expansion [1/K]
thermal_diffusivity_coefficient     = 1.0e-6        # value for the coefficient of thermal diffusivity [m^2/s]
gravity_acceleration                = 10.0          # value for the gravity acceleration [m/s^2]
density_mantle                      = 3300.         # value for the mantle reference density [kg/m^3]
external_heat                       = 0.0e-12       # ok
heat_capacity                       = 1250.         # value for the heat capacity [J/K]

non_linear_method                   = on            # set if non linear method is used for the momentum equation [on/off]
adiabatic_component                 = on            # set if adiabatic heating/cooling is active [on/off]
radiogenic_component                = on            # set if radiogenic heating is active [on/off]

# Velocity boundary conditions
top_normal_velocity                 = fixed         # set the normal velocity on the top side of the model to be fixed or free [fixed/free] 
top_tangential_velocity             = free          # set the tangential velocity on the top side of the model to be fixed or free [fixed/free]
bot_normal_velocity                 = fixed         # set the normal velocity on the bottom side of the model to be fixed or free [fixed/free]
bot_tangential_velocity             = free          # set the tangential velocity on the bot side of the model to be fixed or free [fixed/free]
left_normal_velocity                = fixed         # set the normal velocity on the left side of the model to be fixed or free [fixed/free]
left_tangential_velocity            = fixed         # set the tangential velocity on the left side of the model to be fixed or free [fixed/free]
right_normal_velocity               = fixed         # set the normal velocity on the right side of the model to be fixed or free [fixed/free]
right_tangential_velocity           = fixed         # set the tangential velocity on the right side of the model to be fixed or free [fixed/free]

surface_velocity                    = 0.0e-2        # ok

# Temperature boundary conditions
top_temperature                     = fixed         # set temperature on the top side of the model to be fixed or free [fixed/free]
bot_temperature                     = fixed         # set temperature on the bottom side of the model to be fixed or free [fixed/free]
left_temperature                    = fixed         # set temperature on the left side of the model to be fixed or free [fixed/free]
right_temperature                   = fixed         # set temperature on the right side of the model to be fixed or free [fixed/free]

rheology_model                      = 9             # flag n. of a pre-defined rheology model to use during simulation

# Create the parameter file
with open("param.txt", "w") as f:
    for line in params.split("\n"):
        line = line.strip()
        if len(line):
            f.write(" ".join(line.split()) + "\n")

Run the model

In this example, mandyoc use the following flags:

-seed 0,2
-strain_seed 0.0,1.0

the first flag -seed specify which layers will have a specific initial strain. These layers are the number 0 and 2 for the asthenosphere and the mantle weak seed in the middle of the lithospheric mantle.

The second flag -strain_seed specify the magnitude of the initial strain, with respect to each layer indicated in the previous flag. In these case, the asthenosphere will have no initial strain and the mantle seed will have a initial constant strain of 1.0. The other layers will present a random perturbation for the value of the initial strain, and the maximum magnitude of this initial strain is defined by the parameter random_initial_strain in parameters file.

You can run the model as:

mpirun -n NUMBER_OF_CORES mandyoc -seed 0,2 -strain_seed 0.0,1.0

Where the NUMBER_OF_CORES value must be set by the user.

This process took around 2 days to finish the temporal evolution.


After Mandyoc finished running the model, different files are obtained with the temporal evolution of the model as it was specified in the parameter file. In this case, we are only going to use the density, cumulative strain, temperature, velocity, viscosity and density files.

Plot the results

Now, we will plot the Mandyoc output file for only 5 different steps times. These files are into output/ directory.

Determine the initial and final step to make the plots:

step_initial = 0
step_final = 5000

d_step = 1000

The dark and light orange colors represent the upper and lower crust, respectively, while the dark and light green colors represent lithospheric and sublithospheric mantle, respectively. Shades of gray indicate the magnitude of the cumulative strain. The blue and orange bars indicate the width of the extended continental crust in both conjugate margins.

for cont in range(step_initial, step_final + d_step, d_step):  
    # Read the files:
    # Time
    filename = "output/time_" + str(cont) + ".txt"
    time = load_time(filename)
    # Density
    filename = "output/density_" + str(cont) + ".txt"
    rho = load_data(filename, nx, nz)    
    # Strain
    filename = "output/strain_" + str(cont) + ".txt"
    strain = load_data(filename, nx, nz)    
    strain[rho < 200] = 0
    strain_log = np.log10(strain)
    # Temperature
    filename = "output/temperature_" + str(cont) + ".txt"
    temp = load_data(filename, nx, nz)    
    # Velocity
    filename = "output/velocity_" + str(cont) + ".txt"
    veloc_x, veloc_z = load_velocity(filename, nx, nz)
    # Density
    filename = "output/viscosity_" + str(cont) + ".txt"
    visco = load_data(filename, nx, nz)
    visco_log = np.log10(visco)
    visco_log[visco_log == -np.inf] = 0
    # Create the plots:
    fig, axs = plt.subplots(nrows=3, ncols=1, figsize=(20, 10), sharex=True)
    axs[0].set_title("Time = %.1lf Myr\n\n" % (time[0]/1.0e6))
    # Plot density and strain_log
    # Density
    img_rho = axs[0].contourf(
        (zz + h_sair)/1e3,
        levels=[100, 2700, 2800, 3354, 3378],
        label="Density [$kg/m^3$]"
    # Strain_log
    strain_aux = strain_log.copy()
    strain_aux[strain_aux == -np.Inf] = np.nan
    strain_aux[strain_aux < -0.5] = np.nan
    strain_aux[strain_aux > 0.9] = np.nan
    img_strain = axs[0].scatter(
        (zz + h_sair)/1e3,
    axs[0].set_ylabel("Depth [km]")

    # Plot temperature and velocity
    # Temperature
    img_temp = axs[1].pcolormesh(
        label="Temperature [°C]"
    # Velocity
        xx[::20, ::20]/1e3, 
        (zz[::20, ::20] + h_sair)/1e3,
        (veloc_x[::20, ::20]/1e3).T,
        (veloc_z[::20, ::20]).T/1e3,
    axs[1].set_ylabel("Depth [km]")
    # Plot viscosity
    img_visco = axs[2].pcolormesh(
        label="Viscosity [Pa s]"
    axs[2].set_ylabel("Depth [km]")
    axs[2].set_xlabel("x [km]")  
<Figure size 1440x720 with 6 Axes><Figure size 1440x720 with 6 Axes><Figure size 1440x720 with 6 Axes><Figure size 1440x720 with 6 Axes><Figure size 1440x720 with 6 Axes><Figure size 1440x720 with 6 Axes>

Interesting links

  1. Silva, R. M., & Sacek, V. (2022). Influence of Surface Processes on Postrift Faulting During Divergent Margins Evolution. Tectonics, 41(2). 10.1029/2021tc006808
  2. Salazar-Mora, C. A., & Sacek, V. (2021). Lateral flow of thick continental lithospheric mantle during tectonic quiescence. Journal of Geodynamics, 145, 101830. 10.1016/j.jog.2021.101830
  3. Sacek, V. (2017). Post-rift influence of small-scale convection on the landscape evolution at divergent continental margins. Earth and Planetary Science Letters, 459, 48–57. 10.1016/j.epsl.2016.11.026
  4. Sacek, V., Assunção, J., Pesce, A., & da Silva, R. M. (2022). Mandyoc: A finite element code to simulate thermochemical convection in parallel. Journal of Open Source Software, 7(71), 4070. 10.21105/joss.04070