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)

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

First, let's load our modules:

In [1]:
import gplately
import numpy as np
import pygplates
import glob, os
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from plate_model_manager import PlateModelManager

Define get_plate_tectonic_stats and obtain plate tectonic stats for the Muller et al. 2019 model.

In [2]:
pm_manager = PlateModelManager()
muller2019_model = pm_manager.get_model("Muller2019", data_dir="plate-model-repo")

rotation_model = muller2019_model.get_rotation_model()
topology_features = muller2019_model.get_topologies()
static_polygons = muller2019_model.get_static_polygons()

coastlines = muller2019_model.get_layer('Coastlines')
continents = muller2019_model.get_layer('ContinentalPolygons')
COBs =  muller2019_model.get_layer('COBs')

model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)

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):

  1. longitude of sampled ridge point
  2. latitude of sampled ridge point
  3. spreading velocity magnitude (in cm/yr)
  4. 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):

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

In [3]:
# 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)
    
    # 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
    )
    
    return data

Let's obtain Muller2019 plate tectonic stats from 250Ma to present day.

In [4]:
# Create a time array
min_time = 0 
max_time = 249 #time 250 is no good
time_array = np.arange(max_time,min_time-1, -1)

# Initialise data ndarray
data_array = np.zeros((9,time_array.size))

# Get plate tectonic stats for each reconstruction time
for t, time in enumerate(time_array):
    data = get_plate_tectonic_stats(model, time)
    if data:
        data_array[:,t] = data
    gplately.tools.update_progress((max_time-time)/(max_time-min_time))
print("Done calculating stats!")

Progress: [####################] 100.0%
Done calculating stats!

Save plate tectonic stats to a CSV file¶

In [5]:
headers = ['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)']

# Turn data array into a pandas dataframe.
muller19_df = pd.DataFrame(np.column_stack(data_array), columns = headers)

# Save dataframe to CSV
muller19_df.to_csv('./NotebookFiles/Muller2019-PlateTectonicStats.csv', encoding='utf-8', index=False)

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.

In [6]:
# Smooth data with a Gaussian filter
from scipy.ndimage import gaussian_filter

# Get stats from the data frame
reconstruction_times = muller19_df['Time (Ma)'].to_list()
ridge_vels_mean = muller19_df['Mean global ridge velocities (cm/yr)'].to_list()
ridge_vels_std = muller19_df['Global ridge velocity standard deviation (cm/yr)'].to_list()
subd_vels_mean = muller19_df['Mean global subduction velocities (cm/yr)'].to_list()
subd_vels_std = muller19_df['Global subduction velocity standard deviation (cm/yr)'].to_list()

# Use a Gaussian filter
ridge_vels_mean_smoothed = gaussian_filter(ridge_vels_mean, sigma=2)
subd_vels_mean_smoothed = gaussian_filter(subd_vels_mean, sigma=2)

# Plotting functions
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)
Out[6]:
<matplotlib.legend.Legend at 0x1ee48a41010>
No description has been provided for this image

Calculate plate tectonic stats of other plate models¶

We repeat the workflow above to calculate the plate tectonic stats of two additional plate reconstruction models:

  1. Muller et al. 2016
  2. Merdith 2021
In [7]:
# MULLER ET AL 2016
# Create a time array
min_time = 0 
max_time = 230
muller16_time_array = np.arange(max_time,min_time-1, -1)

# Initialise data ndarray
muller16_data = np.zeros((9,muller16_time_array.size))

# Create the PlateReconstruction object with Muller et al 2016
gdownload = gplately.download.DataServer("Muller2016")
rotation_model, topology_features, static_polygons = gdownload.get_plate_reconstruction_files()
muller2016 = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)

# Get plate tectonic stats for Muller et al 2016
for t, time in enumerate(muller16_time_array):
    muller16_data[:,t] = get_plate_tectonic_stats(muller2016, time)
    #print("Calculated stats for {} Ma!".format(time))
    gplately.tools.update_progress((max_time-time)/(max_time-min_time))
print("Done calculating stats!")

Progress: [####################] 100.0%
Done calculating stats!
In [8]:
# Turn the Muller et al 2016 plate tectoncic stats array into a pandas dataframe.
muller16_df = pd.DataFrame(np.column_stack(muller16_data), columns = headers)

# Save dataframe to CSV
muller16_df.to_csv("./NotebookFiles/Muller2016-PlateTectonicStats.csv", encoding='utf-8', index=False)
In [9]:
# MERDITH 2021
# Create a time array
min_time = 0 
max_time = 250
merdith21_time_array = np.arange(max_time,min_time-1, -1)

# Initialise data ndarray
merdith21_data = np.zeros((9,merdith21_time_array.size))

# Create the PlateReconstruction object with Merdith 2021
gdownload = gplately.download.DataServer("Merdith2021")
rotation_model, topology_features, static_polygons = gdownload.get_plate_reconstruction_files()
merdith2021 = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)

# Get plate tectonic stats for Merdith 2021
for t, time in enumerate(merdith21_time_array):
    merdith21_data[:,t] = get_plate_tectonic_stats(merdith2021, time)
    #print("Calculated stats for {} Ma!".format(time))
    gplately.tools.update_progress((max_time-time)/(max_time-min_time))
print("Done calculating stats!")

Progress: [####################] 100.0%
Done calculating stats!
In [10]:
# Turn the Merdith 2021 plate tectoncic stats array into a pandas dataframe.
merdith21_df = pd.DataFrame(np.column_stack(merdith21_data), columns = headers)

# Save dataframe to CSV
merdith21_df.to_csv("./NotebookFiles/Merdith2021-PlateTectonicStats.csv", encoding='utf-8', index=False)

Visualising plate tectonic stats¶

Read and plot data from the plate tectonic stats dataframes we just created for the Muller et al. 2016 and Merdith 2021 plate models using the pandas library.

Global mid-ocean ridge and subduction zone lengths¶

In [11]:
fig, axes = plt.subplots(2, 1, figsize=(8,6), dpi=200)
fig.subplots_adjust(hspace=0.05)
n_axes = len(axes)

# Access subduction zone and ridge lengths from pandas dataframes
muller19_trench_lengths = muller19_df['Global subduction zone lengths (km)'].to_list()
muller16_trench_lengths = muller16_df['Global subduction zone lengths (km)'].to_list()
merdith21_trench_lengths = merdith21_df['Global subduction zone lengths (km)'].to_list()
muller19_ridge_lengths = muller19_df['Global mid-ocean ridge lengths (km)'].to_list()
muller16_ridge_lengths = muller16_df['Global mid-ocean ridge lengths (km)'].to_list()
merdith21_ridge_lengths = merdith21_df['Global mid-ocean ridge lengths (km)'].to_list()


# Prepare different time arrays for each model as Muller 2016 stats start at 230 Ma, while 
# the others start at 250 Ma
muller19_time_array = muller19_df['Time (Ma)'].to_list()
muller16_time_array = muller16_df['Time (Ma)'].to_list()
muller21_time_array = merdith21_df['Time (Ma)'].to_list()


# Subduction zone length plot
axes[0].set_ylabel('Subduction zone\nlengths (km)')
axes[0].set_title("(a)", x=-0.1, y=1, fontsize="12")
axes[0].plot(muller19_time_array, gaussian_filter(muller19_trench_lengths, sigma=1), c='k', linestyle=":", label="Muller et al. 2019")
axes[0].plot(muller16_time_array, gaussian_filter(muller16_trench_lengths,sigma=1), c='k', alpha = 0.7, label='Muller et al. 2016')
axes[0].plot(muller21_time_array, gaussian_filter(merdith21_trench_lengths,sigma=1), c='k', alpha = 0.3, label="Merdith 2021")
axes[0].tick_params(direction="in", length=5, top=True, right=True)
axes[0].tick_params(direction="in", which="minor", length=2.5, top=True, right=True)
axes[0].minorticks_on()

# Mid-ocean ridge length plot
axes[1].set_ylabel('Mid-ocean ridge\nlengths (km)')
axes[1].set_title("(b)", x=-0.1, y=1, fontsize="12")
axes[1].plot(muller19_time_array, gaussian_filter(muller19_ridge_lengths,sigma=1), c='k', linestyle=":", label="Muller et al. 2019")
axes[1].plot(muller16_time_array, gaussian_filter(muller16_ridge_lengths,sigma=1), c='k', alpha = 0.7, label="Muller et al. 2016")
axes[1].plot(muller21_time_array, gaussian_filter(merdith21_ridge_lengths,sigma=1), c='k', alpha = 0.3, label="Merdith 2021")
axes[1].legend(bbox_to_anchor=(0.9, -0.2), ncol = 3, frameon=False)
axes[1].tick_params(direction="in", length=5, top=True, right=True)
axes[1].tick_params(direction="in", which="minor", length=2.5, top=True, right=True)
axes[1].minorticks_on()

# Place tick labels only on the bottom plot
for i, ax in enumerate(axes):
    ax.set_xlim(max_time, min_time)
    if i < n_axes-1:
        ax.set_xticklabels([])
    else:
        ax.set_xlabel("Age (Ma)")
No description has been provided for this image

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

In [12]:
fig, axes = plt.subplots(2, 1, figsize=(8,6), dpi=200)
fig.subplots_adjust(hspace=0.05)
n_axes = len(axes)

# Access subduction convergence and ridge spreading rates (along with their standard deviation) from the dataframes
# Turn lists into numpy arrays.
muller19_subduction_rate = np.array(muller19_df['Mean global ridge velocities (cm/yr)'].to_list())
muller16_subduction_rate = np.array(muller16_df['Mean global ridge velocities (cm/yr)'].to_list())
merdith21_subduction_rate = np.array(merdith21_df['Mean global ridge velocities (cm/yr)'].to_list())
muller19_subduction_rate_sd = np.array(muller19_df['Global ridge velocity standard deviation (cm/yr)'].to_list())
muller16_subduction_rate_sd = np.array(muller16_df['Global ridge velocity standard deviation (cm/yr)'].to_list())
merdith21_subduction_rate_sd = np.array(merdith21_df['Global ridge velocity standard deviation (cm/yr)'].to_list())

muller19_spreading_rate = np.array(muller19_df['Mean global subduction velocities (cm/yr)'].to_list())
muller16_spreading_rate = np.array(muller16_df['Mean global subduction velocities (cm/yr)'].to_list())
merdith21_spreading_rate = np.array(merdith21_df['Mean global subduction velocities (cm/yr)'].to_list())
muller19_spreading_rate_sd = np.array(muller19_df['Global subduction velocity standard deviation (cm/yr)'].to_list())
muller16_spreading_rate_sd = np.array(muller16_df['Global subduction velocity standard deviation (cm/yr)'].to_list())
merdith21_spreading_rate_sd = np.array(merdith21_df['Global subduction velocity standard deviation (cm/yr)'].to_list())

# Subduction convergence rates with standard deviation
axes[0].set_ylabel('Ridge \n spreading \n rate (cm/yr)', rotation="horizontal", 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])
# Muller et al. 2019
axes[0].plot(muller19_time_array, gaussian_filter(muller19_subduction_rate, sigma=1), c='k', linestyle=":", label="Muller et al. 2019")
axes[0].fill_between(muller19_time_array, 
                 gaussian_filter(muller19_subduction_rate-muller19_subduction_rate_sd, sigma=2), 
                 gaussian_filter(muller19_subduction_rate+muller19_subduction_rate_sd, sigma=2),
                 edgecolor='k', color='k', alpha=0.05)
# Muller et al. 2016
axes[0].plot(muller16_time_array, gaussian_filter(muller16_subduction_rate,sigma=1), c='k', alpha = 0.7, label='Muller et al. 2016')
axes[0].fill_between(muller16_time_array, 
                 gaussian_filter(muller16_subduction_rate-muller16_subduction_rate_sd, sigma=2), 
                 gaussian_filter(muller16_subduction_rate+muller16_subduction_rate_sd, sigma=2),
                 edgecolor='k', color='k', alpha=0.05)
# Merdith 2021
axes[0].plot(merdith21_time_array, gaussian_filter(merdith21_subduction_rate,sigma=1), c='k', alpha = 0.3, label="Merdith et al. 2021")
axes[0].fill_between(merdith21_time_array, 
                 gaussian_filter(merdith21_subduction_rate-merdith21_subduction_rate_sd, sigma=2), 
                 gaussian_filter(merdith21_subduction_rate+merdith21_subduction_rate_sd, sigma=2),
                 edgecolor='k', color='k', alpha=0.05)
# Tick labels
axes[0].tick_params(direction="in", length=5, top=True, right=True)
axes[0].tick_params(direction="in", length=5, top=True, right=True)
axes[0].tick_params(direction="in", which="minor", length=2.5, top=True, right=True)
axes[0].minorticks_on()


# Mid-ocean ridge spreading rates with standard deviation
axes[1].set_ylabel('Subduction \n convergence \n rate (cm/yr)', rotation="horizontal", 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])
axes[1].plot(muller19_time_array, gaussian_filter(muller19_spreading_rate,sigma=1), c='k', linestyle=":", label="Muller et al. 2019")
axes[1].fill_between(muller19_time_array, 
                 gaussian_filter(muller19_spreading_rate-muller19_spreading_rate_sd, sigma=2), 
                 gaussian_filter(muller19_spreading_rate+muller19_spreading_rate_sd, sigma=2),
                 edgecolor='k', color='k', alpha=0.05)

axes[1].plot(muller16_time_array, gaussian_filter(muller16_spreading_rate,sigma=1), c='k', alpha = 0.7, label="Muller et al. 2016")
axes[1].fill_between(muller16_time_array, 
                 gaussian_filter(muller16_spreading_rate-muller16_spreading_rate_sd, sigma=2), 
                 gaussian_filter(muller16_spreading_rate+muller16_spreading_rate_sd, sigma=2),
                 edgecolor='k', color='k', alpha=0.05)

axes[1].plot(merdith21_time_array, gaussian_filter(merdith21_spreading_rate,sigma=1), c='k', alpha = 0.3, label="Merdith et al. 2021")
axes[1].fill_between(merdith21_time_array, 
                 gaussian_filter(merdith21_spreading_rate-merdith21_spreading_rate_sd, sigma=2), 
                 gaussian_filter(merdith21_spreading_rate+merdith21_spreading_rate_sd, sigma=2),
                 edgecolor='k', color='k', alpha=0.05)

axes[1].legend(bbox_to_anchor=(0.9, -0.25), ncol = 3, frameon=False)
axes[1].tick_params(direction="in", length=5, top=True, right=True)
axes[1].tick_params(axis="x", which="major", direction="inout", length=5)
axes[1].tick_params(direction="in", which="minor", length=2.5, top=True, right=True)
axes[1].minorticks_on()

# Place tick labels on the bottom plot only.
for i, ax in enumerate(axes):
    ax.set_xlim(max_time, min_time)
    if i < n_axes-1:
        ax.set_xticklabels([])
    else:
        ax.set_xlabel("Age (Ma)", labelpad=10)
No description has been provided for this image

Global rates of crustal production and destruction¶

In [13]:
fig, axes = plt.subplots(2, 1, figsize=(8,6), dpi=200)
fig.subplots_adjust(hspace=0.05)
n_axes = len(axes)


# Access crustal production and destruction rates from the dataframes
# Turn lists into numpy arrays.
muller19_crustal_production_rate = np.array(muller19_df['Surface area of crust produced at ridges (km^2/yr)'].to_list())
muller16_crustal_production_rate = np.array(muller16_df['Surface area of crust produced at ridges (km^2/yr)'].to_list())
merdith21_crustal_production_rate = np.array(merdith21_df['Surface area of crust produced at ridges (km^2/yr)'].to_list())

muller19_crustal_destruction_rate = np.array(muller19_df['Surface area of crust subducted at trenches (km^2/yr)'].to_list())
muller16_crustal_destruction_rate = np.array(muller16_df['Surface area of crust subducted at trenches (km^2/yr)'].to_list())
merdith21_crustal_destruction_rate = np.array(merdith21_df['Surface area of crust subducted at trenches (km^2/yr)'].to_list())


# Crustal production rates
axes[0].set_ylabel('Crustal \n production \n rate (km$^2$/yr)', rotation="horizontal", 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].plot(muller19_time_array, gaussian_filter(muller19_crustal_production_rate, sigma=1), c='k', linestyle=":", label="Muller et al. 2019")
axes[0].plot(muller16_time_array, gaussian_filter(muller16_crustal_production_rate,sigma=1), c='k', alpha = 0.7, label='Muller et al. 2016')
axes[0].plot(merdith21_time_array, gaussian_filter(merdith21_crustal_production_rate,sigma=1), c='k', alpha = 0.3, label="Merdith et al. 2021")

# Tick labels
axes[0].tick_params(direction="in", length=5, top=True, right=True)
axes[0].tick_params(direction="in", length=5, top=True, right=True)
axes[0].tick_params(direction="in", which="minor", length=2.5, top=True, right=True)
axes[0].minorticks_on()


# Crustal destruction rates
axes[1].set_ylabel('Crustal \n destruction \n rate (km$^2$/yr)', rotation="horizontal", 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].plot(muller19_time_array, gaussian_filter(muller19_crustal_destruction_rate,sigma=1), c='k', linestyle=":", label="Muller et al. 2019")
axes[1].plot(muller16_time_array, gaussian_filter(muller16_crustal_destruction_rate,sigma=1), c='k', alpha = 0.7, label="Muller et al. 2016")
axes[1].plot(merdith21_time_array, gaussian_filter(merdith21_crustal_destruction_rate,sigma=1), c='k', alpha = 0.3, label="Merdith et al. 2021")

axes[1].legend(bbox_to_anchor=(0.9, -0.25), ncol = 3, frameon=False)
axes[1].tick_params(direction="in", length=5, top=True, right=True)
axes[1].tick_params(axis="x", which="major", direction="inout", length=5)
axes[1].tick_params(direction="in", which="minor", length=2.5, top=True, right=True)
axes[1].minorticks_on()

# Place tick labels on the bottom plot only.
for i, ax in enumerate(axes):
    ax.set_xlim(max_time, min_time)
    if i < n_axes-1:
        ax.set_xticklabels([])
    else:
        ax.set_xlabel("Age (Ma)", labelpad=10)

fig.savefig("crustal_production_destruction.pdf", bbox_inches='tight')
No description has been provided for this image
In [ ]: