7 - Working with Plate Tectonic Stats

In this notebook, we use GPlately to calculate plate tectonic stats like the:

• Total length of all mid-ocean ridges (km)
• Mean ridge spreading velocity (spreading rate) (cm/yr)
• Mean ridge spreading velocity standard deviation (cm/yr)
• Crustal surface area produced over 1 yr at ridges (km^2/yr)
• Total length of all subduction zones (km)
• Mean subduction velocity (convergence rate) (cm/yr)
• Mean subduction velocity standard deviation (cm/yr)
• Crustal surface area subducted by trenches over 1 yr (km^2/yr)
• Length of all transform boundaries (km)
• Mean global trench velocities (cm/yr)
• Global trench velocity standard deviation (cm/yr)

for a given plate model (a GPlately PlateReconstruction object) and a reconstruction_time.

import gplately
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter
from plate_model_manager import PlateModelManager

data_dir = Path("WorkflowData/07-Working-with-Plate-Tectonic-Stats")
output_dir = data_dir / "output"
output_dir.mkdir(parents=True, exist_ok=True)

Define get_plate_tectonic_stats below, then use it to calculate plate tectonic stats for all four plate reconstruction models used in this notebook:

1. Müller et al. 2019
2. Müller et al. 2016
3. Merdith et al. 2021
4. Zahirovic et al. 2022

Calculate plate tectonic stats

Calculate statistics at uniformly spaced points along mid-ocean ridges and subduction zones.

For mid-ocean ridges we use gplately.PlateReconstruction.tessellate_mid_ocean_ridges. This uniformly samples along mid-ocean ridge plate boundaries and, for each sample point, returns the following statistics (by default):

0. longitude of sampled ridge point
1. latitude of sampled ridge point
2. spreading velocity magnitude (in cm/yr)
3. length of arc segment (in degrees) that current point is on

Note: This is the default output. You can optionally output extra data such as the spreading obliquity.

Note: The transform segments of mid-ocean ridges are ignored (by default). In other words, sample points (and hence statistics) are not generated along the transform segments. And transform segments are determined by how much the spreading direction deviates from the segment normal (with the default threshold deviation being an empirically determined value that works well for the common global models).

For subduction zones we use gplately.PlateReconstruction.tessellate_subduction_zones. This uniformly samples along subduction zone plate boundaries and, for each sample point, returns the following statistics (by default):

0. longitude of sampled trench point
1. latitude of sampled trench point
2. subducting convergence (relative to trench) velocity magnitude (in cm/yr)
3. subducting convergence velocity obliquity angle in degrees (angle between trench normal vector and convergence velocity vector)
4. trench absolute (relative to anchor plate) velocity magnitude (in cm/yr)
5. trench absolute velocity obliquity angle in degrees (angle between trench normal vector and trench absolute velocity vector)
6. length of arc segment (in degrees) that current point is on
7. trench normal (in subduction direction, ie, towards overriding plate) azimuth angle (clockwise starting at North, ie, 0 to 360 degrees) at current point
8. subducting plate ID
9. trench plate ID

Note: This is the default output. You can optionally output extra data such as the distance to the nearest edge of the trench.

To get more accurate values for the global crustal production and destruction rates we use gplately.PlateReconstruction.crustal_production_destruction_rate. Unlike the above two functions, this function only considers convergence and divergence along plate boundaries. In other words, it does not care whether a plate boundary is labelled as a mid-ocean ridge or a subduction zone (and it includes plate boundaries that are not mid-ocean ridges or subduction zones). Essentially if a location on a plate boundary is diverging then it produces crust, and if it's diverging then it destroys crust. This function, crustal_production_destruction_rate, then sums all diverging points to get the global crustal production rate, and sums all converging points to get the global crustal destruction rate.

Note: An alternative way to calculate the crustal production rate is to explicitly sum the ridge spreading velocity along mid-ocean ridges (obtained using tessellate_mid_ocean_ridges). And similarly, an alternative way to calculate the crustal destruction rate is to explicitly sum the orthogonal subducting convergence velocity along subduction zones (obtained using tessellate_subduction_zones). This is shown in the code below. However, crustal_production_destruction_rate is more accurate as noted above - and it is also shown in the code below (where it overwrites the less accurate values).

# Function to obtain plate tectonic stats at a specific reconstruction time
def get_plate_tectonic_stats(model, reconstruction_time):

    # -----------------------   MID-OCEAN RIDGES   ------------------------------
    # Calculate mid ocean ridge stats with GPlately
    ridge_data = model.tessellate_mid_ocean_ridges(reconstruction_time)

    if ridge_data is None:
        return

    # Ignore data of ridge segments with negative velocities
    ridge_data = ridge_data[ridge_data[:, 2] >= 0]

    # Latitudes and longitudes of points along ridge segments
    ridge_lon = ridge_data[:, 0]
    ridge_lat = ridge_data[:, 1]

    # Global mid-ocean ridge length at this reconstruction time (km)
    ridge_len = (
        sum(
            np.radians(ridge_data[:, 3])
            * gplately.tools.geocentric_radius(ridge_data[:, 1])
        )
        * 1e-3
    )

    # Mean ridge spreading velocity + its standard deviation (in cm/year)
    ridge_vel = ridge_data[:, 2]  # spreading velocities of ridge segments in cm/yr
    ridge_vel_mean = np.mean(
        ridge_vel
    )  # mean global spreading velocity amongst all ridge segments in cm/yr
    ridge_vel_std = np.std(ridge_vel)  # standard deviation

    # Ridge surface area (km^2/yr)
    # Convert ridge velocities from cm/yr to m/yr; ridge lengths are already in km
    ridge_surface_area = ridge_vel_mean * 1e-5 * ridge_len

    # -----------------------   SUBDUCTION ZONES   ------------------------------
    # Calculate subduction convergence stats with GPlately
    subduction_data = model.tessellate_subduction_zones(reconstruction_time)

    # Latitudes and longitudes of points along trench segments
    subduction_lon = subduction_data[:, 0]
    subduction_lat = subduction_data[:, 1]

    # calculate geocentric earth radius from latitude
    earth_radius = gplately.tools.geocentric_radius(subduction_lat)

    # Global subduction zone length at this reconstruction time (km)
    subduction_len = np.sum(np.radians(subduction_data[:, 6]) * earth_radius) * 1e-3

    # Ensure convergence velocities are positive
    subduction_data[:, 2] = np.clip(subduction_data[:, 2], 0.0, 1e99)

    # Multiply convergence velocities by the cosine of the subduction obliquity angle to get
    # subduction velocities (cm/yr)
    subduction_vel = np.fabs(subduction_data[:, 2]) * np.cos(
        np.radians(subduction_data[:, 3])
    )

    # Global mean subduction velocity and stdev of all trench segments (cm/year)
    subd_vel_mean, subd_vel_std = np.mean(subduction_vel), np.std(subduction_vel)

    # Absolute velocity(cm/year) of the trench’s motion in the orthogonal direction towards the overriding plate.
    # Negative if moving towards the overriding plate (Trench advance)
    # Positive if moving away from the overriding plate (Trench retreat)
    # The purpose of "-np.fabs()" is to allow "np.cos()" to produce the correct sign of value(negative or positive).
    trench_absolute_vel = -np.fabs(subduction_data[:, 4]) * np.cos(
        np.radians(subduction_data[:, 5])
    )

    # Global mean and standard deviation of trench velocity(cm/year) of all trench segments
    trench_abs_vel_mean, trench_abs_vel_std = np.mean(trench_absolute_vel), np.std(
        trench_absolute_vel
    )

    # Area subducted by trenches over 1 yr (km^2/yr)
    # Convert subduction velocities from cm/yr to km/yr; trench lengths are already in km.
    subd_surface_area = subd_vel_mean * 1e-5 * subduction_len

    # Use gplately.PlateReconstruction.crustal_production_destruction_rate() to get a more
    # accurate value for global crustal production/destruction rate.
    (
        total_crustal_production_rate_km_2_per_yr,
        total_crustal_destruction_rate_km_2_per_yr,
    ) = model.crustal_production_destruction_rate(time)
    ridge_surface_area = total_crustal_production_rate_km_2_per_yr
    subd_surface_area = total_crustal_destruction_rate_km_2_per_yr

    # Return a set of boundary stats in a tuple
    data = (
        reconstruction_time,
        ridge_len,
        ridge_vel_mean,
        ridge_vel_std,
        ridge_surface_area,
        subduction_len,
        subd_vel_mean,
        subd_vel_std,
        subd_surface_area,
        trench_abs_vel_mean,
        trench_abs_vel_std,
    )

    return data

Compute stats for each plate model, and save to CSV

Each model supports a different maximum reconstruction time (Müller et al. 2019, for example, only extends back to 249 Ma, not 250 Ma). We loop over all four models, use get_plate_tectonic_stats to build a full time series of statistics for each, and save each one to its own CSV file.

# Plate reconstruction models used in this notebook, and the maximum time (Ma) each one supports.
MODEL_CONFIGS = {
    "muller19": {
        "pmm_name": "Muller2019",
        "label": "Muller et al. 2019",
        "max_time": 249,
    },
    "muller16": {
        "pmm_name": "Muller2016",
        "label": "Muller et al. 2016",
        "max_time": 230,
    },
    "merdith21": {
        "pmm_name": "Merdith2021",
        "label": "Merdith et al. 2021",
        "max_time": 250,
    },
    "zahirovic22": {
        "pmm_name": "Zahirovic2022",
        "label": "Zahirovic et al. 2022",
        "max_time": 250,
    },
}

# Columns match the order of the tuple returned by get_plate_tectonic_stats()
STAT_COLUMNS = [
    "Time (Ma)",
    "Global mid-ocean ridge lengths (km)",
    "Mean global ridge velocities (cm/yr)",
    "Global ridge velocity standard deviation (cm/yr)",
    "Surface area of crust produced at ridges (km^2/yr)",
    "Global subduction zone lengths (km)",
    "Mean global subduction velocities (cm/yr)",
    "Global subduction velocity standard deviation (cm/yr)",
    "Surface area of crust subducted at trenches (km^2/yr)",
    "Mean global trench velocities (cm/yr)",
    "Global trench velocity standard deviation (cm/yr)",
]

pmm = PlateModelManager()
model_dfs = {}  # key -> DataFrame of plate tectonic stats, one entry per model
done_labels = []

for key, cfg in MODEL_CONFIGS.items():
    pmm_model = pmm.get_model(cfg["pmm_name"], data_dir="plate-model-repo")
    reconstruction = gplately.PlateReconstruction(
        pmm_model.get_rotation_model(),
        pmm_model.get_topologies(),
        pmm_model.get_static_polygons(),
    )

    time_array = np.arange(cfg["max_time"], -1, -1)
    data_array = np.zeros((len(STAT_COLUMNS), time_array.size))
    for t, time in enumerate(time_array):
        stats = get_plate_tectonic_stats(reconstruction, time)
        if stats:
            data_array[:, t] = stats

        frac = (cfg["max_time"] - time) / cfg["max_time"] if cfg["max_time"] else 1.0
        gplately.tools.update_progress(frac)

        status = ", ".join(done_labels) if done_labels else "none yet"
        print(f"Working on: {cfg['label']} | Completed: {status}")

    df = pd.DataFrame(np.column_stack(data_array), columns=STAT_COLUMNS)
    df.to_csv(
        output_dir / f"{cfg['pmm_name']}-PlateTectonicStats.csv",
        encoding="utf-8",
        index=False,
    )
    model_dfs[key] = df
    done_labels.append(cfg["label"])

# Final summary
print("\nFinished computing plate tectonic stats for:")
for label in done_labels:
    print(f"- {label}")

Progress: [####################] 100.0%
Working on: Zahirovic et al. 2022 | Completed: Muller et al. 2019, Muller et al. 2016, Merdith et al. 2021

Finished computing plate tectonic stats for:
- Muller et al. 2019
- Muller et al. 2016
- Merdith et al. 2021
- Zahirovic et al. 2022

Plotting plate tectonic stats

As an example, let's plot the ridge spreading rate and trench convergence rate for Müller et al. 2019 as a function of time.

muller19_df = model_dfs["muller19"]

reconstruction_times = muller19_df["Time (Ma)"]
ridge_vels_mean = muller19_df["Mean global ridge velocities (cm/yr)"]
ridge_vels_std = muller19_df["Global ridge velocity standard deviation (cm/yr)"]
subd_vels_mean = muller19_df["Mean global subduction velocities (cm/yr)"]
subd_vels_std = muller19_df["Global subduction velocity standard deviation (cm/yr)"]

ridge_vels_mean_smoothed = gaussian_filter(ridge_vels_mean, sigma=2)
subd_vels_mean_smoothed = gaussian_filter(subd_vels_mean, sigma=2)

fig = plt.figure(figsize=(8, 4), dpi=200)
ax1 = fig.add_subplot(
    111,
    xlim=(250, 0),
    ylim=(0, 18),
    xlabel="Age (Ma)",
    ylabel="Rate (cm/yr)",
    title="Subduction zone convergence and ridge spreading rates (cm/yr):\nMüller et al 2019",
)

ax1.plot(reconstruction_times, ridge_vels_mean_smoothed, label="Ridge spreading rate")
ax1.fill_between(
    reconstruction_times,
    gaussian_filter(ridge_vels_mean_smoothed - ridge_vels_std, sigma=2),
    gaussian_filter(ridge_vels_mean_smoothed + ridge_vels_std, sigma=2),
    edgecolor="k",
    color="C0",
    alpha=0.2,
)

ax1.plot(
    reconstruction_times, subd_vels_mean_smoothed, label="Subduction convergence rate"
)
ax1.fill_between(
    reconstruction_times,
    gaussian_filter(subd_vels_mean_smoothed - subd_vels_std, sigma=2),
    gaussian_filter(subd_vels_mean_smoothed + subd_vels_std, sigma=2),
    edgecolor="k",
    color="C1",
    alpha=0.2,
)
plt.legend(loc="upper right", frameon=False)
plt.show()

Visualising plate tectonic stats

Read and plot data from the plate tectonic stats dataframes we just created for all four plate models, using the pandas library.

Global mid-ocean ridge and subduction zone lengths

# Consistent per-model plotting style, reused across all comparison plots below.
MODEL_STYLE = {
    "muller19": dict(color="k", linestyle=":", alpha=1.0, label="Muller et al. 2019"),
    "muller16": dict(color="k", linestyle="-", alpha=0.7, label="Muller et al. 2016"),
    "merdith21": dict(color="k", linestyle="-", alpha=0.3, label="Merdith et al. 2021"),
    "zahirovic22": dict(
        color="k", linestyle="-.", alpha=0.5, label="Zahirovic et al. 2022"
    ),
}


def plot_model_stat(ax, column, sigma=1, sd_column=None):
    """Plot `column` from every model's dataframe on `ax`, using a consistent style per model.

    If `sd_column` is given, also draws a shaded +/- 1 standard deviation band per model.
    """
    for key, df in model_dfs.items():
        time = df["Time (Ma)"]
        values = df[column]
        ax.plot(time, gaussian_filter(values, sigma=sigma), **MODEL_STYLE[key])
        if sd_column is not None:
            sd = df[sd_column]
            ax.fill_between(
                time,
                gaussian_filter(values - sd, sigma=2),
                gaussian_filter(values + sd, sigma=2),
                edgecolor="k",
                color="k",
                alpha=0.05,
            )


def style_axis(ax, bottom=False):
    """Apply the tick styling used consistently across the comparison plots below."""
    ax.tick_params(direction="in", length=5, top=True, right=True)
    if bottom:
        ax.tick_params(axis="x", which="major", direction="inout", length=5)
    ax.tick_params(direction="in", which="minor", length=2.5, top=True, right=True)
    ax.minorticks_on()


def finalize_time_axes(axes, max_time, labelpad=None):
    """Share the x-axis across stacked panels: hide tick labels except on the bottom panel."""
    for i, ax in enumerate(axes):
        ax.set_xlim(max_time, 0)
        if i < len(axes) - 1:
            ax.set_xticklabels([])
        else:
            ax.set_xlabel("Age (Ma)", labelpad=labelpad)

fig, axes = plt.subplots(2, 1, figsize=(8, 6), dpi=200)
fig.subplots_adjust(hspace=0.05)

axes[0].set_ylabel("Subduction zone\nlengths (km)")
axes[0].set_title("(a)", x=-0.1, y=1, fontsize="12")
plot_model_stat(axes[0], "Global subduction zone lengths (km)")
style_axis(axes[0])

axes[1].set_ylabel("Mid-ocean ridge\nlengths (km)")
axes[1].set_title("(b)", x=-0.1, y=1, fontsize="12")
plot_model_stat(axes[1], "Global mid-ocean ridge lengths (km)")
axes[1].legend(bbox_to_anchor=(0.9, -0.2), ncol=3, frameon=False)
style_axis(axes[1], bottom=True)

max_time = max(df["Time (Ma)"].max() for df in model_dfs.values())
finalize_time_axes(axes, max_time)
plt.show()

Global mid-ocean ridge spreading rates and subduction zone convergence rates

fig, axes = plt.subplots(3, 1, figsize=(8, 9), dpi=200)
fig.subplots_adjust(hspace=0.05)

axes[0].set_ylabel("Ridge spreading \n rate (cm/yr)", fontsize="12")
axes[0].yaxis.set_label_coords(-0.125, 0.45)
axes[0].set_title("(a)", x=-0.2, y=0.9, fontsize="12")
axes[0].set_ylim([0, 18])
plot_model_stat(
    axes[0],
    "Mean global ridge velocities (cm/yr)",
    sd_column="Global ridge velocity standard deviation (cm/yr)",
)
style_axis(axes[0])

axes[1].set_ylabel("Subduction convergence \n rate (cm/yr)", fontsize="12")
axes[1].yaxis.set_label_coords(-0.125, 0.45)
axes[1].set_title("(b)", x=-0.2, y=0.9, fontsize="12")
axes[1].set_ylim([0, 18])
plot_model_stat(
    axes[1],
    "Mean global subduction velocities (cm/yr)",
    sd_column="Global subduction velocity standard deviation (cm/yr)",
)
style_axis(axes[1])

axes[2].set_ylabel("Trench Velocities (cm/yr)", fontsize="12")
axes[2].yaxis.set_label_coords(-0.125, 0.45)
axes[2].set_title("(c)", x=-0.2, y=0.9, fontsize="12")
axes[2].set_ylim([-9.5, 13.5])
plot_model_stat(
    axes[2],
    "Mean global trench velocities (cm/yr)",
    sd_column="Global trench velocity standard deviation (cm/yr)",
)
axes[2].legend(bbox_to_anchor=(0.9, -0.25), ncol=3, frameon=False)
style_axis(axes[2], bottom=True)

max_time = max(df["Time (Ma)"].max() for df in model_dfs.values())
finalize_time_axes(axes, max_time, labelpad=10)
plt.show()

Global rates of crustal production and destruction

fig, axes = plt.subplots(2, 1, figsize=(8, 6), dpi=200)
fig.subplots_adjust(hspace=0.05)

axes[0].set_ylabel("Crustal production \n rate (km$^2$/yr)", fontsize="12")
axes[0].yaxis.set_label_coords(-0.125, 0.45)
axes[0].set_title("(a)", x=-0.2, y=0.9, fontsize="12")
plot_model_stat(axes[0], "Surface area of crust produced at ridges (km^2/yr)")
style_axis(axes[0])

axes[1].set_ylabel("Crustal destruction \n rate (km$^2$/yr)", fontsize="12")
axes[1].yaxis.set_label_coords(-0.125, 0.45)
axes[1].set_title("(b)", x=-0.2, y=0.9, fontsize="12")
plot_model_stat(axes[1], "Surface area of crust subducted at trenches (km^2/yr)")
axes[1].legend(bbox_to_anchor=(0.9, -0.25), ncol=3, frameon=False)
style_axis(axes[1], bottom=True)

max_time = max(df["Time (Ma)"].max() for df in model_dfs.values())
finalize_time_axes(axes, max_time, labelpad=10)

fig.savefig(output_dir / "crustal_production_destruction.pdf", bbox_inches="tight")
plt.show()