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

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.


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(Müller et al., , 2008a(Müller et al., , 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).

Methods
We divide the computational approach for generating seafloor age grids from kinematic plate models into three modules: preprocessing, main algorithm and a postprocessing (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.

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 timestep. The collection of all the plate velocity fields constitute a global surface velocity field V(t) (Fig. 1A) at each time-step t ∈ [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).  (Matthews et al., 2016).

Plate Velocity [cm/yr]
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 X i , where i = 0, 1, .., N − 1, and N is the number of mesh nodes. Thus, at some time t, a velocity field 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 X i are static, while the properties (velocities 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).
We extract plate boundaries labelled mid-ocean ridge and sample these circle arc segments at a user-defined resolution given by the parameter ∆ 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 (∆ R ∼ 20 − 50 km). The collection of points with the desired spacing define midocean ridge locations, and is denoted 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 S k (t) for k = 0, 1, .., N S (t) − 1, with the resolution given by ∆ S (values of ∼ 15 − 30 km are suggested).
We developed python-scripts that set up and extract the necessary input variables (X i , V i (t), P i (t), R j (t) and S k (t)) to our algorithm from full-plate tectonic reconstructions. These can be found in: http://doi.org/10.5281/zenodo.3383154.

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 x n (t), their age τ 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 ∈ [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

Pre-processing
Post-processing These constitute the input variables: 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.
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 ∆ R , and the age τ 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 X i to the tracer positions 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 X i . Therefore V i (t) must be interpolated to the tracer positions 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 x n (t) are denoted v n (t).

4) Move tracers
To obtain the next positions of the tracers (Fig. 3 B-C) at the time-step t + ∆t, we use the forward-Euler method: x n (t + ∆t) = x n (t) + ∆tv n (t). The default timestep ∆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 ∼ 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 + ∆t) are interpolated to the new tracer positions x n (t + ∆t), obtaining p n (t + ∆t).

6) Check for subduction
After having moved the tracers from x n (t) to x n (t + ∆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 + ∆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 closedpolygon 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 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 ∆t to their current age, obtaining τ n (t + ∆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).

Post-Processing
The output from steps 1-8) is randomly distributed tracer positions x n (t) associated with ages τ 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.

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 x n (t = 0) and τ 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.

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)

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 ±1 Myr of the ages from our age grids, and about 80 % are within ±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 ∼ 45 % and ∼ 92 % respectively (Fig. S2C), and ∼ 52 % and ∼ 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. Karlsen et al. 2019a), which should preclude a constant agearea 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 timedependent age-area distributions of the seafloor as generated by our algorithm with those published in Müller et al. (2016). 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 algorithmgenerated 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.

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 agedepth 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(Worsley et al., , 1986Nance 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.

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 sub-duction 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.

A) B) C)
Paleozoic Mesozoic Cenozoic Pangea 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.

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;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.