Running the model

To run the accompanying .prm files simulating rifting in a continental setting, execute the following lines with ASPECT:

%%capture
! aspect-release prm_files/rifting.prm

In the above .prm file, rifting is initiated due to the preexisting weakness in the plastic_strain compositional field (described in the previous notebook.

%%capture
! aspect-release prm_files/rifting_with_upwelling.prm

On the other hand, the above .prm file simulates rifting in presence of an elevated thermal anomaly.

In this above cells, it is assumed that the ASPECT directory is installed system-wide or that you are using the HubZero ASPECT installation. If not, modify the ASPECT executable to the location where it is installed.

ASPECT can be configured in either “release” mode or “debug” mode. Release mode is more optimized than debug mode, and will usually run much faster but lacks debugging symbols and more descriptive error messages. If you are making changes to ASPECT or altering parameter files, it might be a good idea to try running in debug mode in order to resolve any problems. The above code is calling the release version using aspect-release . If you want to try running debug mode, call aspect-debug instead.

Note: The provided .prm files are configured to run on a single node in a reasonable time, but a finer mesh is recommended that can resolve the rift formations. To do this, increase `set Initial global refinement` to a higher value. The plots below use 3 as the global refinement level.

import os
os.environ['VTK_DEFAULT_OPENGL_WINDOW'] = 'vtkEGLRenderWindow' # Suppresses a warning with Pyvista
# Load the relevant libraries
import pyvista as pv
import glob
from IPython.display import HTML

pv.global_theme.allow_empty_mesh = True

def plot_deformation(output_dir, field, gif_file, cmap, model_title, clim):
    '''
    Process a series of VTU files in the given output directory, returning the 'field'
    data for plotting.
    Inputs:
    -----------
    output_dir : str
        Path to the folder containing solution VTU files (expected in solution/solution-*.vtu)
    field : str
        Name of the field we want to look at in the solution
    label : str
        Label for the field being plotted
    cmap  : str
        Colorscale to use for plotting the field
    model_title : str
        Title to be displayed on the plot
    clim : list
        Color limits for the plot (e.g., [0, 1])
    '''

    solution_file_names = sorted(glob.glob(output_dir + 'solution/solution-*.pvtu'))

    mesh = pv.read(solution_file_names[0])

    plotter = pv.Plotter()

    # Create the first frame
    plotter.add_mesh(mesh, scalars=field, show_scalar_bar=False)
    plotter.add_axes()

    plotter.hide_axes()

    # We know that the model is XY plane so the following
    # lines just sets up the camera
    plotter.view_xy()
    plotter.camera.roll = 0
    plotter.camera.parallel_projection = True

    sargs_1 = dict(height=0.5, vertical=True, position_x=0.05, position_y=0.05, n_labels=7)

    plotter.open_gif(gif_file, fps=4)

    # Update function for each frame
    def update(frame):
        frame_int = int(frame)
        new_mesh = pv.read(solution_file_names[frame_int])

        plotter.add_mesh(new_mesh, scalars=field, cmap=cmap, scalar_bar_args=sargs_1, \
                         clim=clim, log_scale=True)

        plotter.add_text(model_title)

    n_frames = len(solution_file_names)

    for i in range(n_frames):
        update(i)
        plotter.write_frame()

    plotter.close()

plot_deformation('rifting-with-upwelling/', 'strain_rate', 'images/rifting_with_upwelling_strain_rate.gif', 'magma', 'Rifting with thermal upwelling', [1e-20, 1e-14])
plot_deformation('rifting/', 'strain_rate', 'images/rifting_strain_rate.gif', 'magma', 'Rifting with distributed weak seeds', [1e-20, 1e-14])

Now, you would have the gifs saved in the current directory this cell visualizes the model outputs side by side for comparison.

In the above model setups, rifts begin to form at the locations of weak seeds in the lithosphere (left) or distributed uniformly in the crust (right) in the absence of pre-existing weaknesses. Notice that the high strain rates in the upwelling regions are associated with low viscosities due to temperature-dependence of the rheology.

Exercises

• How does the model evolve if you change the extension rate? Does the surface topography change?

• Increase the strength of the lithosphere by increasing the activation energy in the rheology. How does this affect the rift evolution?

Appendix

The example file provided takes about 5 minutes on 4 cores. When the model is complete, you should expect the following output on your terminal: