# A Tracer-Based Algorithm for Automatic Generation of Seafloor Age Grids from
Plate Tectonic Reconstructions

Krister S. Karlsen (k.s.karlsen@geo.uio.no), Mathew Domeier,
Carmen Gaina and Clinton P. Conrad

Centre for Earth Evolution and Dynamics, University of Oslo, Norway

Abstract

The age of the ocean floor and its time-dependent age distribution control
fundamental features of the Earth, such as bathymetry, sea level and mantle
heat loss. Recently, the development of increasingly sophisticated
reconstructions of past plate motions has provided models for plate kinematics
and plate boundary evolution back in geological time. These models implicitly
include the information necessary to determine the age of ocean floor that has
since been lost to subduction. However, due to the lack of an automated and
efficient method for generating global seafloor age grids, many tectonic
models, most notably those extending back into the Paleozoic, are published
without an accompanying set of age models for oceanic lithosphere. Here we
present an automatic, tracer-based algorithm that generates seafloor age grids
from global plate tectonic reconstructions with defined plate boundaries. Our
method enables us to produce the first seafloor age models for the Paleozoic’s
lost ocean basins. Estimated changes in sea level based on bathymetry inferred
from our new age grids show good agreement with sea level record estimations
from proxies, providing a possible explanation for the peak in sea level
during the assembly phase of Pangea. This demonstrates how our seafloor age
models can be directly compared with observables from the geologic record that
extend further back in time than the constraints from preserved seafloor.
Thus, our new algorithm may also aid the further development of plate tectonic
reconstructions by strengthening the links between geological observations and
tectonic reconstructions of deeper time.

##  1 Introduction

The discovery of a method to determine the age of the present-day oceanic
crust, using reversals of the Earth’s magnetic field (Vine and Matthews,
1963), gave rise to the recognition that the seafloor is spreading, and
ultimately to the development and broad acceptance of plate tectonics. In the
half-century since the plate tectonic revolution, detailed age models of the
present-day oceanic lithosphere have been constructed, and digital global
oceanic age grids are continuously refined (Müller et al., 1997, 2008a, 2016).
A wealth of information, mainly from marine geophysical data, but also from
geology of continental margins, were used to reconstruct the extent and age
distribution of oceanic lithosphere of the past, including portions that have
been subducted (Müller et al., 2008b). These “paleo-seafloor age grids”
present rich new opportunities for scientific inquiry, as a wide range of
Earth processes can be further interrogated with the use of such age grids.
Example applications include the estimation of paleobathymetry (spatial and
temporal changes in ocean basin depth, which in turn is important for
understanding past ocean currents and their effect on paleoclimate, e.g.,
Straume et al. 2019), sea level change (Müller et al., 2008b), global seafloor
heat flow (Loyd et al., 2007; Crameri et al., 2019), and the subduction volume
flux, which impacts geomagnetic reversals (Hounslow et al., 2018), the thermal
structure of paleo-subduction zones (Maunder et al., 2019), water transport to
the deep mantle (Karlsen et al., 2019a), and the slab pull force on tectonic
plates (Conrad and Lithgow-Bertelloni, 2004; Faccenna et al., 2012). Seafloor
ages for past times are also important as a boundary condition for global
mantle convection models (Gurnis et al., 2012).

Present day age grid models are based on a set of isochrons (lines defined by
equal seafloor ages) constructed using information from magnetic and gravity
data available at various resolutions in most oceanic basins. Ages for
seafloor locations between isochrons are computed based on rotation parameters
that describe the plate motions for various time intervals. The isochrons and
rotation parameters are linked to a specific geomagnetic timescale, and the
choice of timescales will influence the calculated values of spreading
velocities. To ensure a smooth grid of ocean floor ages that maintains sharp
age discontinuities at fracture zones, Müller et al. (1997) designed an
algorithm where they first created a set of densely interpolated isochrons
along plate flow lines, and then used a minimum curvature routine to obtain
age values on a regular grid. This method for reconstructing seafloor age from
present-day seafloor age data and a plate kinematic model is unfortunately
extremely time-consuming and requires significant human input, and
consequently may be subjective or introduce errors. Because of this, seafloor
ages are usually only determined after a plate reconstruction model has been
finalized; they have not previously been computed "on the fly" from the plate
kinematic model itself.

The mid-ocean ridges constitute the locus of seafloor generation through time,
while plate kinematics define the seafloor’s subsequent journey until its
destruction at a subduction zone. Thus, global plate tectonic reconstructions
that define the motions of the plates and the locations and types of plate
boundaries (for a detailed description of this type of model see Gurnis et al.
2012) also implicitly define the age and structure of oceanic basins. Global
plate tectonic models with dynamic plate boundaries have been constructed back
into the mid-Paleozoic (410 Ma; Domeier and Torsvik 2014; Matthews et al.
2016), but published paleo-seafloor age grids are only available globally for
the last 250 Ma (Müller et al., 2019). This timing discrepancy has, partly,
occurred because the reconstruction of global paleo-seafloor age grids from a
plate tectonic reconstruction presents a tedious and labor intensive task in
the absence of an automatic method. Moreover, these reconstructions are
subject to continuous changes as new geological information becomes available.
It follows that an automatic and more efficient method for seafloor age
determination is needed to allow the use of plate tectonic reconstructions to
better decipher Earth’s processes and dynamics through time. An automatic
method also allows detecting inconsistencies in kinematic models, which would
help to improve them. In this study we present such an automatic method for
generating seafloor age grids, and introduce a specific implementation of it
as an open-source Python code called _Tracer Tectonics_ (or _TracTec_)
(Karlsen et al., 2019b).

##  2 Methods

We divide the computational approach for generating seafloor age grids from
kinematic plate models into three modules: preprocessing, main algorithm and a
post-processing (Fig. 2). This division saves time by avoiding expensive data
processing along the way, so that the main algorithm steps can remain as
streamlined and efficient as possible.

###  2.1 Preprocessing

Global, full-plate tectonic reconstructions present a set of dynamically
closed plate polygons that evolve over time from a set of dynamically evolving
plate boundaries (Gurnis et al., 2012), which are derived from a wealth of
geological and geophysical data. Each plate is associated with a _Plate ID_
and its motion can be defined by rotation about an Euler pole, or an
equivalent surface velocity field, at each time-step. The collection of all
the plate velocity fields constitute a _global surface velocity field_
$\mathbf{V}(t)$ (Fig. 1A) at each time-step $t\in[0,T]$, where $t=0$ defines
the present-day and $t=T$ is the earliest point in time for which the plate
reconstruction is defined. The dynamically evolving plate boundaries are
labelled according to the type of plate boundary that they represent; here we
are interested in only those that are either a mid-ocean ridge or a subduction
zone (Fig. 1B).

Figure 1: Present day global surface velocity field (A) and locations of mid-
ocean ridges and subduction zones (B) derived from global plate tectonic
reconstructions (Matthews et al., 2016).

As the time-dependent spatial distribution of mid-ocean ridges and subduction
zones, together with plate kinematics, dictates the age of the ocean floor, we
need to extract these properties from a given full-plate model and output them
in a convenient format. To accomplish this, we use _Pygplates_ ; a Python-
based scripting interface to GPlates (Boyden et al., 2011) that allows for
easy automation of such tasks. To describe the plate kinematics, we define a
global mesh with node coordinates given by $\mathbf{X}_{i}$, where
$i=0,1,..,N-1$, and $N$ is the number of mesh nodes. Thus, at some time $t$, a
velocity field $\mathbf{V}_{i}(t)$ and a set of Plate IDs $P_{i}(t)$ can be
assigned to each node of the mesh by interpolation from the full-plate model.
Additionally, to be able to distinguish between oceanic and continental
regions, we also extract the locations of continents through time and
interpolate them to our globally defined mesh, such that each mesh node is
associated with a binary region value $C_{i}(t)$ through time. $C_{i}(t)=1$
indicates that this mesh node is associated with a continent, while
$C_{i}(t)=0$ represents an ocean. Note that the mesh node coordinates
$\mathbf{X}_{i}$ are static, while the properties (velocities
$\mathbf{V}_{i}(t)$, Plate IDs $P_{i}(t)$ and region type $C_{i}(t)$) at these
nodes change with time. There are no requirement on the type of mesh, for
instance it could be a regular latitude-longitude mesh, or a CitcomS mesh
(Zhong et al., 2000).

Figure 2: Flow chart of TracTec (Karlsen et al., 2019b) showing an overview of
the algorithm and the steps used to generate seafloor age grids from plate
tectonic reconstructions.

We extract plate boundaries labelled mid-ocean ridge and sample these circle
arc segments at a user-defined resolution given by the parameter $\Delta_{R}$,
defined the distance (km) between each point (or vertex). Because this
resolution parameter also controls the number of tracers that are added at the
mid-ocean ridges to track the evolution of the seafloor age (see next
section), we recommend a relatively high resolution ($\Delta_{R}\sim 20-50$
km). The collection of points with the desired spacing define mid-ocean ridge
locations, and is denoted $\mathbf{R}_{j}(t)$, where $j=0,1,..,N_{R}(t)-1$ and
$N_{R}(t)$ is the number of ridge vertices at time $t$. Subduction zone plate
boundary segments are extracted and sampled in the same manner, and the
vertices that define their locations are given by $\mathbf{S}_{k}(t)$ for
$k=0,1,..,N_{S}(t)-1$, with the resolution given by $\Delta_{S}$ (values of
$\sim 15-30$ km are suggested).

We developed python-scripts that set up and extract the necessary input
variables ($\mathbf{X}_{i}$, $\mathbf{V}_{i}(t)$, $P_{i}(t)$,
$\mathbf{R}_{j}(t)$ and $\mathbf{S}_{k}(t)$) to our algorithm from full-plate
tectonic reconstructions. These can be found in:
http://doi.org/10.5281/zenodo.3383154.

###  2.2 Main Algorithm

To simulate and track the journey of oceanic lithosphere through space and
time; from its creation along a ridge until its destruction at a subduction
zone, we use tracers. These are numerical particles on which quantities of
interest are tracked. The essential property we want to study is the age of
the oceanic lithosphere, and this is the primary property we track with the
tracers. A secondary property we track is the Plate ID associated with each
tracer, i.e. the plate to which each tracer currently belongs. We track Plate
IDs to determine when a tracer crosses a plate boundary (recognized by a
change in its associated Plate ID), which in turn is useful to detect
subduction of tracers. At a given time $t$, the number of tracers is
$N_{T}(t)$, their positions are $\mathbf{x}_{n}(t)$, their age $\tau_{n}(t)$,
and their associated Plate IDs are $p_{n}(t)$, for $n=0,1,..,N_{T}(t)-1$.

We break the algorithm into individual sub-steps (see below, 1-8) that are
completed at every time-step $t\in[0,T]$, before moving on to the next:

_1) Add tracers_
At the beginning of each time-step we add tracers on each side of the mid-
ocean ridges, with an offset distance of 30 km away from the ridge axis to
each side. This is to ensure that tracers get assigned the correct plate
velocity from the global velocity field. To demonstrate the importance of this
strategy we provide an example: the lowermost ridge-axis in Figure 3 A-C lines
up parallel to the mesh nodes. If tracers were added exactly (zero offset
distance) along the ridge-axis, all would be assigned a velocity to the left,
and a sparse region would open up as they were moved. The spacing between the
added tracers (on each side of the ridge) is given by $\Delta_{R}$, and the
age $\tau_{n}(t)$ of these tracers is set to one (this is equivalent of
assuming an average spearing rate of 3 cm/yr for the first one Myr).

_2) Interpolate Plate ID to tracers_
Plate IDs $P_{i}(t)$ are interpolated from the mesh nodes $\mathbf{X}_{i}$ to
the tracer positions $\mathbf{x}_{n}(t)$ using the N-dimensional nearest
neighbor algorithm from the Python library Scipy (Jones et al., 2001),
obtaining $p_{n}(t)$.

_3) Interpolate velocity to tracers_
The tracers that track the oceanic lithosphere’s journey must move according
to the plate velocities known only at the mesh nodes $\mathbf{X}_{i}$.
Therefore $\mathbf{V}_{i}(t)$ must be interpolated to the tracer positions
$\mathbf{x}_{n}(t)$ (Fig. 3B). We employ the same nearest neighbor algorithm
for the velocities, because this conserves the discontinuous nature of
tectonic plate velocities and avoids unphysical smoothing across plate
boundaries, while providing superior efficiency on a sphere. The interpolated
velocities at $\mathbf{x}_{n}(t)$ are denoted $\mathbf{v}_{n}(t)$.

_4) Move tracers_
To obtain the next positions of the tracers (Fig. 3 B-C) at the time-step
$t+\Delta t$, we use the forward-Euler method: $\mathbf{x}_{n}(t+\Delta
t)=\mathbf{x}_{n}(t)+\Delta t\mathbf{v}_{n}(t)$. The default time-step $\Delta
t$ is set to 1 Myr, which is appropriately small in the context of average
plate speeds. Notably, although our algorithm is perfectly suited to operate
at any arbitrarily shorter (or larger) time-steps, the temporal resolution of
most full-plate models is $\sim$ 1 Myr. Thus, use of shorter time-steps can
result in erroneous behavior owing to shortcomings in the temporal resolution
of the input model.

_5) Update Plate IDs_
Plate IDs $P_{i}(t+\Delta t)$ are interpolated to the new tracer positions
$\mathbf{x}_{n}(t+\Delta t)$, obtaining $p_{n}(t+\Delta t)$.

_6) Check for subduction_
After having moved the tracers from $\mathbf{x}_{n}(t)$ to
$\mathbf{x}_{n}(t+\Delta t)$, we check if any tracers have been subducted
(Fig. 3 H-I). This is accomplished by first determining which tracers have
changed Plate ID by comparing $p_{n}(t)$ against $p_{n}(t+\Delta t)$. Given
the sub-set of tracers that have changed Plate ID, we determine if any are
within distance $R_{min}$ of a subduction zone. The default value of $R_{min}$
is set to 100 km, which we have found to be appropriate. Any tracers that
fulfill both of these criteria are considered subducted and are thus
terminated at this time-step.

By default, the subducted tracers are stored in auxiliary files that can be
used for studying various properties of the subduction history, for example to
estimate the average age of subducting plates through time, or to consider the
time-integrated history of subduction.

_7) Check for tracer-continent interaction_
To ensure that tracers don’t end up on continents (which shouldn’t happen in
closed-polygon models, but is not impossible in nature (e.g., ophiolites) and
can occur due to inexorable flaws in full-plate models), we delete all tracers
whose closest mesh node $\mathbf{X}_{i}$ is associated with a continent
$C_{i}(t)$. From the perspective of computational efficiency, it is also
preferred to maintain only the tracers that are strictly needed.

_8) Update tracer ages_
Before moving on to the next time-step, which practically means returning to
step 1) of the algorithm, we update the age of the tracers that are left after
the two filtering steps 6-7), by simply by adding $\Delta t$ to their current
age, obtaining $\tau_{n}(t+\Delta t)$.

Note that as pygplates (and the GPlates files associated with the input plate
model) are used only in preprocessing step, the age grid algorithm (steps 1-8)
could be applied to output from a global mantle convection model, in the case
the user had a way of tracking the distribution of plates and their associated
boundaries and velocities at the surface (e.g. Mallard et al. 2017).

###  2.3 Post-Processing

The output from steps 1-8) is randomly distributed tracer positions
$\mathbf{x}_{n}(t)$ associated with ages $\tau_{n}(t)$. However, for most
applications, it is convenient to express the seafloor ages on a regular grid,
rather than at random points. This calls for a post-processing step that
interpolates seafloor ages from tracer positions onto a regular grid. There
are countless ways of doing this, ranging from simple nearest neighbor
algorithms, to weighted means, to splines etc. Depending on the application,
smoothing of resulting agegrids may be preferred or required. In our online
repository we provide an example of a simple post-processing script that uses
GMT’s linear interpolation algorithm (Wessel et al., 2013) to obtain seafloor
ages on a regular grid.

The Python scripts to generate age grids based on the steps described above
using Tracer Tectonics is provided in the Zenodo repository (Karlsen et al.,
2019b). A summary and a general overview of the algorithm are shown in Figure
2.

###  2.4 Initial Condition

Running the algorithm without any initial condition, as in the example of
Figure 3 D-G, we see that it takes some tens of millions of years before the
ocean basins are entirely covered by tracers. Technically, this is the time it
takes for the predicted seafloor ages to result solely from the plate
kinematics of the input model. Alternatively, one could apply an initial
condition that incorporates some educated guess of the seafloor ages for the
initial time step. This is straightforward to include in the framework
presented (steps 1-8 of the main algorithm), by simply initializing
$\mathbf{x}_{n}(t=0)$ and $\tau_{n}(t=0)$. As with any time-dependent model,
one should be aware of the assumptions implicit to the chosen initial
condition, and its effects on the model output. For this algorithm (and in
contrast to many geodynamic models), it is straightforward to track the effect
of an applied initial condition. This can be achieved by simply tracking the
fraction of initial tracers through time (Fig. S1); at some point, all the
initial tracers will be eliminated, from which point the output will no longer
depend on the initial condition.

##  3 Discussion

###  3.1 Validation and Benchmarking

For computational methods to effectively contribute to advances in the field
of geoscience (as well as in the natural sciences in general), output from
algorithms and simulations should be validated against relevant geological
model constraints. This is particularly true when introducing new modeling
approaches and software, making the practice of code benchmarking an essential
and necessary part of software development. To validate our algorithm and its
implementation, we compare the seafloor age grids generated by our algorithm
against published present-day seafloor age models (Fig. 4), direct point
observations of present-day seafloor ages (Fig. S2) and against the time-
dependent age-area distribution of the oceanic lithosphere for the last 230
Myr (Fig. 5) based on M16 (Müller et al., 2016).

Figure 3: Figure 3: Schematics illustrating how tracers are added along mid-
ocean ridges (A), moved (B-C), and checked for subduction (H-I). The repeated
application of these steps over time leads to an ocean basin gradually filled
with tracers that track the age of the ocean floor (D-G). In this example the
model of Matthews et al. (2016) is used. Note that (A-C, H-I) are not to
scale, thus density of tracers and mesh nodes shown here are not
representative.

We used Matthews et al. (2016) as the input plate model, with the initial
condition shown in Figure 6A, to generate the seafloor age grids used to
evaluate the performance of our algorithm in this section.

To generate the seafloor age grids used to evaluate the performance of our
algorithm in this section, we employed Matthews et al. (2016) as the input
plate model with the initial condition shown in Figure 6A. Notably, however,
the effect of an initial condition at 400 Ma is very small already by 300 Ma,
and zero at the present-day (Fig. S1). The initial condition and the age grids
are available from our online repository (Karlsen et al., 2019b).

Our algorithm reproduces the present-day seafloor well (Fig. 4 A-C), as can be
seen from the maps comparing the resulting age grids to those of M16 (Fig. 4
E-G). The characteristic triangular present-day seafloor age-area distribution
(Sclater et al. 1980; Cogné and Humler 2004), which shows the fraction of the
ocean floor that falls within a certain age range, is also well reproduced. On
the regional scale, some minor differences can be seen, for example there are
three narrow bands of artificially young seafloor branching from the Mid-
Atlantic Ridge near the Caribbean and the south-west of Iberia. These are
merely consequences of the underlying plate model (Matthews et al., 2016), for
which these plate boundaries are erroneously designated as mid-ocean ridges,
either at present (Fig. 1), or in the recent past. This demonstrates that a
combined work-flow for developing plate tectonic reconstructions, which
includes analysis of the seafloor ages, can reveal flaws and inconsistencies
in the model. In the case of the aforementioned errors in the mid-Atlantic,
simply re-labeling the offending boundaries as transform features and re-
running the age grid algorithm addresses the issue.

To further evaluate the performance of our algorithm, we compare the generated
present-day age grid against ages inferred from magnetic reversal picks (Fig.
S2). This global dataset provides by far the most comprehensive, direct
sources of oceanic lithosphere ages, and is available from the Global Seafloor
Fabric and Magnetic Lineation (GSFML) database (Seton et al., 2014). We see
that about 30 % of the 101418 pick ages are within $\pm 1$ Myr of the ages
from our age grids, and about 80 % are within $\pm 5$ Myr (Fig. S2B). These
numbers are slightly higher for models of the present-day ocean floor that are
based on more labor-intensive methods (described in Section 1) such as Muller
et al. (2008b), for which they are $\sim$ 45 % and $\sim$ 92 % respectively
(Fig. S2C), and $\sim$ 52 % and $\sim$ 93 % for M16 (Fig. S2D).

As pointed out in several studies (e.g. Becker et al. 2009; Coltice et al.
2012; Sim et al. 2016), the triangular age-area distribution of the present-
day ocean floor is unlikely to be a constant feature through Earth’s history.
In particular, large fluctuations are predicted to have occurred in the rates
of seafloor spreading and global subduction, (e.g. Domeier et al. 2018;
Karlsen et al. 2019a), which should preclude a constant age-area distribution
of the seafloor. Moreover, the age-area distribution of the seafloor is an
important feature of our planet because it exhibits a first-order control on
e.g. bathymetry, sea level, and planetary cooling through regulation of
surface heat flow (i.e. loss of mantle heat). Therefore, as a third benchmark,
we compare the time-dependent age-area distributions of the seafloor as
generated by our algorithm with those published in Müller et al. (2016).

Figure 4: Maps comparing present-day seafloor ages generated by our algorithm
(A-C) to those of M16 (Müller et al., 2016) model (E-G), and their
corresponding age-area distributions (D and H). The input plate model used
here was (Matthews et al., 2016).

We observe a clear time-dependence in the age-area distributions computed from
our age grids, and this time-dependence is nearly identical to that predicted
by M16 (Fig. 5). Given the broad observed trend toward relatively younger
seafloor ages between 140-80 Ma, we anticipate that some of these seafloor age
variations are related to the widespread development of new ridge systems
during the Cretaceous. The initiation of some of these new ridge systems
(including the southern mid-Atlantic) was associated with late-stage Pangea
breakup, whereas other new ridges appeared in the Paleo-Pacific basin,
probably related to the emplacement of large igneous provinces (Torsvik et
al., 2019). The rapid development of these new ridges, producing juvenile
oceanic lithosphere, must have been synchronously balanced with increased
subduction that would have accelerated the destruction of relatively old
seafloor; and these two processes combined thus explains the tendency toward
relatively younger seafloor during the Cretaceous (Fig. 5).

In summary, our algorithm reproduces detailed reconstructions of seafloor ages
like those of M16. We observe only minor regional differences between the
algorithm-generated present-day seafloor age grid and the M16 model. These
differences mainly occur in regions where Matthews et al. (2016) inferred
plate boundary locations that deviate from the interpreted isochrons of M16.
This merely demonstrates that our algorithm is an automatic and fast way to
detect inconsistencies between data and models. Global features like the age-
area distribution, which is important for many geodynamic applications, are
robustly reproduced over time. Thus, as plate tectonic reconstructions
improve, so will the reliability and regional resolution of the age-grids
produced by our algorithm.

###  3.2 An Example Application: Paleo Sea Level

In this section we will demonstrate how our new algorithm and paleo-age grids
(computed as described in Section 3.1) can be used to study tectonic
mechanisms for sea level variations during the last 400 Myr. We would like to
stress that the age-grids generated by our algorithm are bound to the input
plate model, and as is the case for all plate models operating in times
earlier than 200 Ma, the kinematics of all oceanic plates are necessarily
synthetic (because no in situ oceanic lithosphere older than 200 Ma has
survived to the present-day). This naturally implies that the construction of
age-grids for earlier times is much more uncertain, and interpretation of them
should be done with caution and care; here we only use these synthetic age
grids to consider some global, first-order trends for the sake of
demonstration.

From the generated seafloor age grids we compute bathymetry by applying the
age-depth relation of (Crosby and McKenzie, 2009). Next, we use these
bathymetry grids (e.g., Fig. 6B) to compute the change in average ocean basin
depth relative to the present-day. Finally, we compare how these changes in
average ocean depth would affect sea level, and relate them to the sea level
history (Fig. 6C) as inferred from the sedimentary record (Hallam, 1992; Haq
and Al-Qahtani, 2005; Haq and Schutter, 2008). Although many other processes
affect sea level fluctuations on tectonic time scales, the age-area
distribution of the seafloor (through thermal subsidence of the ageing oceanic
lithosphere) exhibits the first-order control, and is the only process that
shows a direct correlation with the sea level record (Müller et al., 2008b;
Conrad, 2013; Karlsen et al., 2019a). Therefore, an automatic method for
generating seafloor age grids (from which ocean depth can be computed) enables
sea level to be used as a first-order order deep-time constraint on plate
tectonic reconstructions.

Figure 5: Comparison of the age-area distribution of the seafloor through time
for our algorithm created age grids (top) and the (Müller et al., 2016) (M16)
model (bottom).

Our prediction of sea level fluctuations caused by ocean depth changes
inferred from our 400 Myr reconstruction of age grids shows a first-order
agreement with established sea level records (Hallam, 1992; Haq and Al-
Qahtani, 2005; Haq and Schutter, 2008). Such a correlation has been noted
previously into the Mesozoic (e.g., Müller et al. 2008b; Conrad 2013), but our
algorithm applied to the plate reconstructions of Matthews et al. 2016 shows
that this correlation may extend back to the late Paleozoic. Our age models
predict a clear peak in sea level during the late stages of Pangea assembly
(Fig. 6C), in agreement with the early predictions of sea-level change based
on conceptual models of a supercontinental cycle (Worsley et al., 1985, 1986;
Nance et al., 1986). Moreover, the consistency between predicted and observed
sea level indicates that the tectonic rates (seafloor spreading and its
counterpart, seafloor subduction) in the underlying plate tectonic model might
be reasonable for this period of time.

###  3.3 Limitations

The uncertainties in the generated seafloor age grids are directly linked to,
and controlled by, the underlying plate tectonic model. Thus, the generated
age-grids are dependent on global plate velocities and the locations of mid-
ocean ridges and subduction zones, for which the assignment of formal errors
is either impractical or impossible. For these reasons, we expect uncertainty
in the age grids to be proportional to those of the underlying plate tectonic
reconstructions, with negligible additional uncertainty introduced through the
rather straightforward computational steps of our algorithm (Fig. 2). A final
word of caution is that the generated age grids will never be better than the
input plate model. We thus advise users to familiarize themselves with the
uncertainties tied to the underlying plate tectonic reconstructions.

Figure 6: Initial condition (A) used to generate seafloor age grids from the
plate reconstruction of Matthews et al. (2016) from 400 Ma to the present, and
associated bathymetry (B) at 300 Ma, computed by applying Crosby and McKenzie
(2009). From the bathymetry models, we compute changes in mean ocean depth and
isostatically compensate them (Pitman III, 1978) to compare sea level changes
relative to the present day (C, red dashed line) to Phanerozoic sea level
reconstructions (C, blue lines) inferred from the sedimentary record (Hallam,
1992; Haq and Al-Qahtani, 2005; Haq and Schutter, 2008). Thick black bar shows
the approximate duration of Pangea.

##  4 Conclusions

The development of plate tectonic reconstructions in the past decade has been
rapid, with the introduction of powerful new tools (Boyden et al., 2011;
Gurnis et al., 2012; Shephard et al., 2017) and the augmentation of global
kinematic datasets (Seton et al., 2014; Gaina and Jakob, 2019) propelling a
concomitant proliferation of new and emergent full-plate models, of both
regional (Shephard et al., 2013; Zahirovic et al., 2014; Gibbons et al., 2015;
Domeier, 2016; Zahirovic et al., 2016; Domeier, 2018; Torsvik et al., 2019)
and global scope (Seton et al., 2012; Domeier and Torsvik, 2014; Müller et
al., 2016; Matthews et al., 2016; Merdith et al., 2017; Müller et al., 2019).
Part of the power of these latest full-plate models is that they implicitly
provide information on plate kinematics across the entire global (or regional)
surface, and so can be used to derive seafloor age grids to test and evaluate
them, as well as to potentially explore other geophysical questions that
relate to seafloor ages. Unfortunately, while the information needed to
compute seafloor ages is implicitly available, the presently-established
workflows to retrieve and process this information are time-consuming,
laborious and not publicly available. This emphasizes the need for a method
that can automatically generate seafloor age grids from full-plate tectonic
reconstructions.

In this study we have presented an algorithm for generating seafloor age grids
that robustly predicts the known present-day ocean floor ages, as well as a
previously published model for the time-dependent age-area distribution. This
new method produced the first set of oceanic lithosphere age grids that
approximate the first order age distribution of the Paleozoic oceans.
Application of our generated seafloor models to estimate past sea level
changes reveals a general agreement with observations from the sedimentary
stratigraphy record (Hallam, 1992; Haq and Al-Qahtani, 2005; Haq and Schutter,
2008), and provides a possible explanation for a peak in sea level during the
assembly phase of Pangea. We hope that our automated algorithm will enable
such comparisons between age-grid predictions and existing geological
constraints to become a routine procedure of the full-plate model development
process. Such an improved workflow should ultimately lead to better and more
self-consistent combined reconstructions of both plate tectonics and seafloor
ages.

## Computer Code Availability

Name of code: Tracer Tectonics (TracTec)
Developer: Krister S. Karlsen
First available: September 1, 2019.
Hardware required: No requirements.
Program language: Python.
Software dependencies: GMT (only post-processing) and the following Python
libraries: pygplates (only preprocessing), scipy, numpy.
Program size: 860 MB (including example input files, benchmark age grids,
etc.)
Source code: https://doi.org/10.5281/zenodo.3383154

## Acknowledgments

This paper has benefited from discussions with members of the CEED Earth
Modelling group. All calculations in the main algorithm were preformed with
Python (Van Rossum and Drake, 2011), while GMT (Wessel et al., 2013) was used
for post-processing. We thank Fabio Crameri for the development of
perceptually uniform and color-vision-deficiency friendly color maps (Crameri,
2018b), ensuring an accurate representation of the underlying data in the
figures (Crameri, 2018a).

This research was funded by the Research Council of Norway’s Centres of
Excellence Project 223272.

## References

  * Becker et al. (2009) Becker, Thorsten W; Conrad, Clinton P; Buffett, Bruce, and Müller, R Dietmar.  Past and present seafloor age distributions and the temporal evolution of plate tectonic heat transport.  _Earth and Planetary Science Letters_ , 278(3-4):233–242, 2009.
  * Boyden et al. (2011) Boyden, James A; Müller, R Dietmar; Gurnis, Michael; Torsvik, Trond H; Clark, James A; Turner, Mark; Ivey-Law, Hamish; Watson, Robin J, and Cannon, John S.  Next-generation plate-tectonic reconstructions using gplates.  2011\.
  * Cogné and Humler (2004) Cogné, Jean-Pascal and Humler, Eric.  Temporal variation of oceanic spreading and crustal production rates during the last 180 my.  _Earth and Planetary Science Letters_ , 227(3-4):427–439, 2004.
  * Coltice et al. (2012) Coltice, Nicolas; Rolf, Tobias; Tackley, Paul J, and Labrosse, Stéphane.  Dynamic causes of the relation between area and age of the ocean floor.  _Science_ , 336(6079):335–338, 2012.
  * Conrad (2013) Conrad, Clinton P.  The solid earth’s influence on sea level.  _GSA Bulletin_ , 125(7-8):1027–1052, 2013.
  * Conrad and Lithgow-Bertelloni (2004) Conrad, Clinton P and Lithgow-Bertelloni, Carolina.  The temporal evolution of plate driving forces: Importance of “slab suction” versus “slab pull” during the cenozoic.  _Journal of Geophysical Research: Solid Earth_ , 109(B10), 2004.
  * Crameri (2018a) Crameri, Fabio.  Geodynamic diagnostics, scientific visualisation and staglab 3.0.  _Geoscientific Model Development_ , 11(6):2541–2562, 2018a.
  * Crameri (2018b) Crameri, Fabio.  Scientific colour-maps.  Zenodo. http://doi.org/10.5281/zenodo.1243862, 2018b.
  * Crameri et al. (2019) Crameri, Fabio; Conrad, Clinton P; Montési, Laurent, and Lithgow-Bertelloni, Carolina R.  The dynamic life of an oceanic plate.  _Tectonophysics_ , 760:107–135, 2019.
  * Crosby and McKenzie (2009) Crosby, AG and McKenzie, D.  An analysis of young ocean depth, gravity and global residual topography.  _Geophysical Journal International_ , 178(3):1198–1219, 2009.
  * Domeier (2016) Domeier, Mathew.  A plate tectonic scenario for the iapetus and rheic oceans.  _Gondwana Research_ , 36:275–295, 2016.
  * Domeier (2018) Domeier, Mathew.  Early paleozoic tectonics of asia: Towards a full-plate model.  _Geoscience Frontiers_ , 9(3):789–862, 2018.
  * Domeier and Torsvik (2014) Domeier, Mathew and Torsvik, Trond H.  Plate tectonics in the late paleozoic.  _Geoscience Frontiers_ , 5(3):303–350, 2014.
  * Domeier et al. (2018) Domeier, Mathew; Magni, Valentina; Hounslow, Mark W, and Torsvik, Trond H.  Episodic zircon age spectra mimic fluctuations in subduction.  _Scientific reports_ , 8(1):17471, 2018.
  * Faccenna et al. (2012) Faccenna, Claudio; Becker, Thorsten W; Lallemand, Serge, and Steinberger, Bernhard.  On the role of slab pull in the cenozoic motion of the pacific plate.  _Geophysical Research Letters_ , 39(3), 2012.
  * Gaina and Jakob (2019) Gaina, Carmen and Jakob, Johannes.  Global eocene tectonic unrest: Possible causes and effects around the north american plate.  _Tectonophysics_ , 760:136–151, 2019.
  * Gibbons et al. (2015) Gibbons, AD; Zahirovic, S; Müller, RD; Whittaker, JM, and Yatheesh, V.  A tectonic model reconciling evidence for the collisions between india, eurasia and intra-oceanic arcs of the central-eastern tethys.  _Gondwana Research_ , 28(2):451–492, 2015.
  * Gurnis et al. (2012) Gurnis, Michael; Turner, Mark; Zahirovic, Sabin; DiCaprio, Lydia; Spasojevic, Sonja; Müller, R Dietmar; Boyden, James; Seton, Maria; Manea, Vlad Constantin, and Bower, Dan J.  Plate tectonic reconstructions with continuously closing plates.  _Computers & Geosciences_, 38(1):35–42, 2012\.
  * Hallam (1992) Hallam, Anthony.  _Phanerozoic sea-level changes_.  Columbia University Press, 1992.
  * Haq and Al-Qahtani (2005) Haq, Bilal U and Al-Qahtani, Abdul Motaleb.  Phanerozoic cycles of sea-level change on the arabian platform.  _GeoArabia_ , 10(2):127–160, 2005.
  * Haq and Schutter (2008) Haq, Bilal U and Schutter, Stephen R.  A chronology of paleozoic sea-level changes.  _Science_ , 322(5898):64–68, 2008.
  * Hounslow et al. (2018) Hounslow, Mark W; Domeier, Mathew, and Biggin, Andrew J.  Subduction flux modulates the geomagnetic polarity reversal rate.  _Tectonophysics_ , 742:34–49, 2018.
  * Jones et al. (2001) Jones, Eric; Oliphant, Travis; Peterson, Pearu, and others, .  SciPy: Open source scientific tools for Python, 2001.  URL http://www.scipy.org/.
  * Karlsen et al. (2019a) Karlsen, Krister S; Conrad, Clinton P, and Magni, Valentina.  Deep water cycling and sea level change since the breakup of pangea.  _Geochemistry, Geophysics, Geosystems_ , 20(6):2919–2935, 2019a.
  * Karlsen et al. (2019b) Karlsen, Krister S; Domeier, Mathew; Gaina, Carmen, and Conrad, Clinton P.  TracerTectonics (version 1.0).  Zenodo, 2019b.  URL http://doi.org/10.5281/zenodo.3383154.
  * Loyd et al. (2007) Loyd, SJ; Becker, TW; Conrad, CP; Lithgow-Bertelloni, C, and Corsetti, FA.  Time variability in cenozoic reconstructions of mantle heat flow: plate tectonic cycles and implications for earth’s thermal evolution.  _Proceedings of the National Academy of Sciences_ , 104(36):14266–14271, 2007.
  * Mallard et al. (2017) Mallard, C; Jacquet, B, and Coltice, N.  Adopt: A tool for automatic detection of tectonic plates at the surface of convection models.  _Geochemistry, Geophysics, Geosystems_ , 18(8):3197–3208, 2017.
  * Matthews et al. (2016) Matthews, Kara J; Maloney, Kayla T; Zahirovic, Sabin; Williams, Simon E; Seton, Maria, and Mueller, R Dietmar.  Global plate boundary evolution and kinematics since the late paleozoic.  _Global and Planetary Change_ , 146:226–250, 2016.
  * Maunder et al. (2019) Maunder, B; van Hunen, J; Bouilhol, P, and Magni, V.  Modeling slab temperature: A reevaluation of the thermal parameter.  _Geochemistry, Geophysics, Geosystems_ , 20(2):673–687, 2019.
  * Merdith et al. (2017) Merdith, Andrew S; Collins, Alan S; Williams, Simon E; Pisarevsky, Sergei; Foden, John D; Archibald, Donnelly B; Blades, Morgan L; Alessio, Brandon L; Armistead, Sheree; Plavsa, Diana, and others, .  A full-plate global reconstruction of the neoproterozoic.  _Gondwana Research_ , 50:84–134, 2017.
  * Müller et al. (1997) Müller, R Dietmar; Roest, Walter R; Royer, Jean-Yves; Gahagan, Lisa M, and Sclater, John G.  Digital isochrons of the world’s ocean floor.  _Journal of Geophysical Research: Solid Earth_ , 102(B2):3211–3214, 1997.
  * Müller et al. (2008a) Müller, R Dietmar; Sdrolias, Maria; Gaina, Carmen, and Roest, Walter R.  Age, spreading rates, and spreading asymmetry of the world’s ocean crust.  _Geochemistry, Geophysics, Geosystems_ , 9(4), 2008a.
  * Müller et al. (2008b) Müller, R Dietmar; Sdrolias, Maria; Gaina, Carmen; Steinberger, Bernhard, and Heine, Christian.  Long-term sea-level fluctuations driven by ocean basin dynamics.  _science_ , 319(5868):1357–1362, 2008b.
  * Müller et al. (2016) Müller, R Dietmar; Seton, Maria; Zahirovic, Sabin; Williams, Simon E; Matthews, Kara J; Wright, Nicky M; Shephard, Grace E; Maloney, Kayla T; Barnett-Moore, Nicholas; Hosseinpour, Maral, and others, .  Ocean basin evolution and global-scale plate reorganization events since pangea breakup.  _Annual Review of Earth and Planetary Sciences_ , 44:107–138, 2016.
  * Müller et al. (2019) Müller, R Dietmar; Zahirovic, Sabin; Williams, Simon E; Cannon, John; Seton, Maria; Bower, Dan J; Tetley, Michael G; Heine, Christian; Le Breton, Eline; Liu, Shaofeng, and others, .  A global plate model including lithospheric deformation along major rifts and orogens since the triassic.  _Tectonics_ , 2019.
  * Nance et al. (1986) Nance, R Damian; Worsley, Thomas R, and Moody, Judith B.  Post-archean biogeochemical cycles and long-term episodicity in tectonic processes.  _Geology_ , 14(6):514–518, 1986.
  * Pitman III (1978) Pitman III, Walter C.  Relationship between eustacy and stratigraphic sequences of passive margins.  _Geological Society of America Bulletin_ , 89(9):1389–1403, 1978.
  * Sclater et al. (1980) Sclater, JjG; Jaupart, Cr, and Galson, D_.  The heat flow through oceanic and continental crust and the heat loss of the earth.  _Reviews of Geophysics_ , 18(1):269–311, 1980\.
  * Seton et al. (2012) Seton, M; Müller, RD; Zahirovic, S; Gaina, C; Torsvik, T; Shephard, G; Talsma, A; Gurnis, M; Turner, M; Maus, S, and others, .  Global continental and ocean basin reconstructions since 200 ma.  _Earth-Science Reviews_ , 113(3-4):212–270, 2012\.
  * Seton et al. (2014) Seton, Maria; Whittaker, Joanne M; Wessel, Paul; Müller, R Dietmar; DeMets, Charles; Merkouriev, Sergey; Cande, Steve; Gaina, Carmen; Eagles, Graeme; Granot, Roi, and others, .  Community infrastructure and repository for marine magnetic identifications.  _Geochemistry, Geophysics, Geosystems_ , 15(4):1629–1641, 2014.
  * Shephard et al. (2013) Shephard, Grace E; Müller, R Dietmar, and Seton, Maria.  The tectonic evolution of the arctic since pangea breakup: Integrating constraints from surface geology and geophysics with mantle structure.  _Earth-Science Reviews_ , 124:148–183, 2013.
  * Shephard et al. (2017) Shephard, Grace Elizabeth; Matthews, Kara J; Hosseini, Kasra, and Domeier, Mathew.  On the consistency of seismically imaged lower mantle slabs.  _Scientific reports_ , 7(1):10976, 2017.
  * Sim et al. (2016) Sim, Shi J; Stegman, Dave R, and Coltice, Nicolas.  Influence of continental growth on mid-ocean ridge depth.  _Geochemistry, Geophysics, Geosystems_ , 17(11):4425–4437, 2016.
  * Straume et al. (2019) Straume, EO; Gaina, Carmen; Medvedev, S; Hochmuth, Katharina; Gohl, Karsten; Whittaker, Joanne M; Abdul Fattah, R; Doornenbal, JC, and Hopper, John R.  Globsed: Updated total sediment thickness in the world’s oceans.  _Geochemistry, Geophysics, Geosystems_ , 20(4):1756–1772, 2019.
  * Torsvik et al. (2019) Torsvik, Trond H; Steinberger, Bernhard; Shephard, Grace E; Doubrovine, Pavel V; Gaina, Carmen; Domeier, Mathew; Conrad, Clinton P, and Sager, William W.  Pacific-panthalassic reconstructions: Overview, errata and the way forward.  _Geochemistry, Geophysics, Geosystems_ , 20(7):3659–3689, 2019.
  * Van Rossum and Drake (2011) Van Rossum, Guido and Drake, Fred L.  _The python language reference manual_.  Network Theory Ltd., 2011.
  * Vine and Matthews (1963) Vine, Fred J and Matthews, Drummond Hoyle.  Magnetic anomalies over oceanic ridges.  _Nature_ , 199(4897):947–949, 1963.
  * Wessel et al. (2013) Wessel, Paul; Smith, Walter HF; Scharroo, Remko; Luis, Joaquim, and Wobbe, Florian.  Generic mapping tools: improved version released.  _Eos, Transactions American Geophysical Union_ , 94(45):409–410, 2013.
  * Worsley et al. (1986) Worsley, Thomas R; Nance, R Damian, and Moody, Judith B.  Tectonic cycles and the history of the earth’s biogeochemical and paleoceanographic record.  _Paleoceanography_ , 1(3):233–263, 1986.
  * Worsley et al. (1985) Worsley, TR; Moody, JB, and Nance, RD.  Proterozoic to recent tectonic tuning of biogeochemical cycles.  _The carbon cycle and atmospheric CO2: natural variations Archean to present_ , 32:561–572, 1985.
  * Zahirovic et al. (2014) Zahirovic, S; Seton, M, and Müller, RD.  The cretaceous and cenozoic tectonic evolution of southeast asia.  _Solid Earth_ , 5(1):227, 2014.
  * Zahirovic et al. (2016) Zahirovic, Sabin; Matthews, Kara J; Flament, Nicolas; Müller, R Dietmar; Hill, Kevin C; Seton, Maria, and Gurnis, Michael.  Tectonic evolution and deep mantle structure of the eastern tethys since the latest jurassic.  _Earth-Science Reviews_ , 162:293–337, 2016.
  * Zhong et al. (2000) Zhong, Shijie; Zuber, Maria T; Moresi, Louis, and Gurnis, Michael.  Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection.  _Journal of Geophysical Research: Solid Earth_ , 105(B5):11063–11082, 2000.