8 - Predicting slab flux

Estimate slab flux of subducting oceanic lithosphere using plate thickness from lithospheric cooling models (Grose, 2012), convergence velocity, and trench-segment length.

Data packages

A plate reconstruction and corresponding seafloor age grids are required to estimate slab flux. These can be downloaded from https://www.earthbyte.org/gplates-2-3-software-and-data-sets/

The workflow below has been tested with the Clennett et al. (2020) and Müller et al. (2019) reconstructions, and should also work with other reconstruction models.

References

• Grose, C. J. (2012). Properties of oceanic lithosphere: Revised plate cooling model predictions. Earth and Planetary Science Letters, 333–334, 250–264. https://doi.org/10.1016/j.epsl.2012.03.037
• Clennett, E. J., Sigloch, K., Mihalynuk, M. G., Seton, M., Henderson, M. A., Hosseini, K., et al. (2020). A Quantitative Tomotectonic Plate Reconstruction of Western North America and the Eastern Pacific Basin. Geochemistry, Geophysics, Geosystems, 21(8), 1–25. https://doi.org/10.1029/2020GC009117
• Müller, R. D., Zahirovic, S., Williams, S. E., Cannon, J., Seton, M., Bower, D. J., et al. (2019). A Global Plate Model Including Lithospheric Deformation Along Major Rifts and Orogens Since the Triassic. Tectonics, 38(6), 1884–1907. https://doi.org/10.1029/2018TC005462

import gplately
import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage import gaussian_filter
from plate_model_manager import PlateModelManager

# This workflow takes long time to finish. For a quick run, set the flag below to True
quick_run = False
#quick_run = True

if not quick_run:
    time_steps = np.arange(0, 231)
else:
    time_steps = np.array([0, 230])

We compare subduction-zone data between two plate reconstructions: Müller et al. (2016) and Müller et al. (2019).

pm_manager = PlateModelManager()

pmm_model_1 = pm_manager.get_model("Muller2019", data_dir="plate-model-repo")
assert pmm_model_1 is not None, "Failed to load the Muller2019 plate model."
rotation_model_1 = pmm_model_1.get_rotation_model()
topology_features_1 = pmm_model_1.get_topologies()
model_1 = gplately.PlateReconstruction(rotation_model_1, topology_features_1)

pmm_model_2 = pm_manager.get_model("Muller2016", data_dir="plate-model-repo")
assert pmm_model_2 is not None, "Failed to load the Muller2016 plate model."
rotation_model_2 = pmm_model_2.get_rotation_model()
topology_features_2 = pmm_model_2.get_topologies()
model_2 = gplately.PlateReconstruction(rotation_model_2, topology_features_2)

Retrieve kinematic data

We extract present-day subduction-zone kinematics to estimate slab flux from subducting plate age, convergence rate, and trench-segment length.

def get_subduction_volume_rate(pmm_model, model, reconstruction_time):
    # Calculate subduction convergence using GPlately.
    subduction_data = model.tessellate_subduction_zones(
        reconstruction_time, output_subducting_absolute_velocity_components=True
    )

    subduction_lon = subduction_data[:, 0]
    subduction_lat = subduction_data[:, 1]
    subduction_angle = subduction_data[:, 3]
    subduction_norm = subduction_data[:, 7]
    subduction_pid_sub = subduction_data[:, 8]
    subduction_pid_over = subduction_data[:, 9]
    subduction_length = np.radians(
        subduction_data[:, 6]
    ) * gplately.tools.geocentric_radius(subduction_data[:, 1])
    subduction_convergence = (
        np.fabs(subduction_data[:, 2])
        * 1e-2
        * np.cos(np.radians(subduction_data[:, 3]))
    )
    subduction_migration = (
        np.fabs(subduction_data[:, 4])
        * 1e-2
        * np.cos(np.radians(subduction_data[:, 5]))
    )
    subduction_plate_vel = subduction_data[:, 10]

    # Remove entries with effectively negative subduction.
    # This occurs when subduction obliquity is greater than 90 degrees.
    subduction_convergence = np.clip(subduction_convergence, 0, 1e99)

    # Sample the AgeGrid for the current timestep.

    # Returns a Raster object.
    graster = gplately.Raster(
        data=pmm_model.get_raster("AgeGrids", reconstruction_time),
        plate_reconstruction=model,
        extent=(-180, 180, -90, 90),
    )
    graster.fill_NaNs(inplace=True)
    age_interp = graster.interpolate(subduction_lon, subduction_lat)

    subduction_age = age_interp
    thickness = gplately.tools.plate_isotherm_depth(age_interp)

    # Calculate subduction volume rate in m^3/yr.
    subduction_vol_rate = (
        thickness * subduction_length * subduction_convergence
    )  # Integrated along total subduction-segment length.
    subduction_vol_rate *= 1e-9  # Convert m^3/yr to km^3/yr.

    mean_plate_thickness = thickness.mean()
    mean_subduction_segment_length = subduction_length.sum()
    mean_subduction_convergence_rate = subduction_convergence.mean()
    total_subduction_volume_rate = subduction_vol_rate.sum()

    return (
        mean_plate_thickness,
        mean_subduction_segment_length,
        mean_subduction_convergence_rate,
        total_subduction_volume_rate,
    )

use_parallel = False

if not use_parallel:
    # Calculate total subduction volume rate (km^3/yr) at each timestep for each model.

    # Müller 2019 model
    thk2019 = np.zeros(time_steps.size)  # mean plate thickness
    len2019 = np.zeros(time_steps.size)  # subduction-zone length
    vel2019 = np.zeros(time_steps.size)  # mean convergence velocity
    vol2019 = np.zeros(time_steps.size)  # subduction flux

    # Müller 2016 model
    thk2016 = np.zeros(time_steps.size)
    len2016 = np.zeros(time_steps.size)
    vel2016 = np.zeros(time_steps.size)
    vol2016 = np.zeros(time_steps.size)


    for t, time in enumerate(time_steps):
        thk2019[t], len2019[t], vel2019[t], vol2019[t] = get_subduction_volume_rate(
            pmm_model_1, model_1, time
        )
        thk2016[t], len2016[t], vel2016[t], vol2016[t] = get_subduction_volume_rate(
            pmm_model_2, model_2, time
        )

        gplately.tools.update_progress(time / time_steps.size)
    gplately.tools.update_progress(1)

Progress: [####################] 100.0%

Parallel processing

GPlately supports parallel execution across multiple CPU cores. Here we use joblib to run get_subduction_volume_rate across many timesteps efficiently.

Note: On Windows, joblib may run slower in parallel than in serial for this workflow. We therefore disable parallel processing on Windows by setting n_jobs=None.

if use_parallel:
    from joblib import Parallel, delayed
    import platform

    # Use serial processing on Windows (parallel can be slower for this workflow).
    if platform.system() == "Windows":
        n_jobs = None
    else:
        n_jobs = -3  # Use all CPUs except two.

    # Use the Loky backend.
    parallel = Parallel(n_jobs=n_jobs, backend="loky", verbose=1)

    muller_2016_data = parallel(
        delayed(get_subduction_volume_rate)(pmm_model_2, model_2, time)
        for time in time_steps
    )
    muller_2019_data = parallel(
        delayed(get_subduction_volume_rate)(pmm_model_1, model_1, time)
        for time in time_steps
    )

    # Unpack arrays.
    thk2016, len2016, vel2016, vol2016 = np.array(muller_2016_data).T
    thk2019, len2019, vel2019, vol2019 = np.array(muller_2019_data).T

Plot slab flux through time

# Plot the slab-flux time series with light Gaussian smoothing.
fig = plt.figure(figsize=(12, 6), dpi=300)

muller2016_volumes_smoothed = gaussian_filter(vol2016, sigma=1)
muller2019_volumes_smoothed = gaussian_filter(vol2019, sigma=1)
plt.plot(
    time_steps,
    muller2016_volumes_smoothed,
    color="k",
    label="Müller et al. (2016)",
)
plt.plot(
    time_steps,
    muller2019_volumes_smoothed,
    linestyle="--",
    alpha=0.5,
    color="k",
    label="Müller et al. (2019)",
)

# Plot settings
plt.title("Total Subduction Flux Through Time")
plt.xlabel("Time (Ma)")
plt.ylabel("Subduction Flux (km$^3$/yr)")
plt.legend(loc="upper center", bbox_to_anchor=(0.5, -0.1), ncol=2)
plt.show()