Pressure‐Driven Poiseuille Flow Inherited From Mesozoic Mantle Circulation Led to the Eocene Separation of Australia and Antarctica

The separation between Australia and Antarctica occurred during the final stages of the break‐up of Pangea. Reconstructions of the rifting of the Australian plate away from Antarctica show fast spreading rates since Mid‐Eocene (45 Ma). These reconstructions can be used to understand and quantify the forces driving the Australia/Antarctica separation, and to test hypotheses on mechanisms that may be of shallow (i.e., lithosphere) or deep (i.e., mantle) origin. Analytical calculations indicate that plate‐boundary forces are highly unlikely to be a plausible candidate to explain such a separation. Thus, we use a recently developed global coupled models of mantle and lithosphere dynamics, here we show that this event, whose kinematics are reproduced in our models within the bounds of the reconstruction uncertainties, owes to a significant degree to the pressure‐driven asthenospheric Poiseuille flow associated with the mantle buoyancy field inherited from viscous circulation history throughout the Mesozoic. On the contrary, in simulations when such a buoyancy field is replaced by another one resulting from a random distribution of mantle temperature–thus not representative of Earth’s mantle circulation history–the rapid northward motion of Australia does not occur. Similarly, suppressing contemporaneous plate‐boundary processes (i.e., subduction of the Pacific ridge at the Aleutians and healing of the India‐Australia ridge) from our models does not have a noticeable effect on the Australia‐Antarctica kinematics. Thus, a pressure‐driven Poiseuille mantle flow must be considered, at least in this example and possible elsewhere, as a main driver of plate tectonics.

Plate motions and their changes through time represent unique constraints to understand the past and present evolution of flow in the Earth's mantle. Past plate motions are reconstructed mainly from observations of the ocean-floor magnetization, which records the imprints of past reversals of Earth's magnetic field. An accurate representation of these observations is essential to study plate tectonics and mantle dynamics (e.g., Müller et al., 1997). This has motivated efforts towards building global models of past plate motions (e.g., Müller et al., 2008Müller et al., , 2016Torsvik et al., 2010, among others), from which plate-motion changes can be inferred. The current observational coverage of the magnetization of ocean floor younger than ∼20 Myr allows reconstructing paleo-kinematic reconstructions at unprecedented temporal resolution of around 1 Myr (e.g., DeMets & Merkouriev, 2019;DeMets et al., 2015). The record of temporal plate-motion changes stands out, revealing variations that occur on short time-scales relative to the time it takes for the largescale structure associated with mantle buoyancy to evolve (Bunge et al., 1998). Some of the reconstructed plate-motion changes have been linked to evolving plate boundary processes and topographic loads (e.g., Copley et al., 2010;Iaffaldano et al., 2006), while other such changes may owe to sublithospheric processes (e.g., Cande & Stegman, 2011).
In this study we focus on the reconstructed history of divergence between the Australian and Antarctica plates. We aim at understanding the role that shallow-seated (i.e., generated within the lithosphere) and deep-seated (i.e., originating from the mantle buoyancy field) processes might have had in driving the separation between the two plates. To this end, we make use of a recently developed numerical model that solves the time-dependent balance of torques within the coupled mantle and lithosphere system. Our numerical model provides explicit predictions of global plate motions (i.e., in the form of Euler vectors) and associated driving/resisting torques as they change through geological time (Stotz et al., 2017(Stotz et al., , 2018. We perform several numerical simulations in order to (i) reproduce the characteristics of the reconstructed Australia/ Antarctica motion, and (ii) determine whether the Australia/Antarctica divergence is driven by torques originating within the lithosphere, or arising from the mantle large-scale buoyancy field.

Reconstructed Australia/Antarctica Relative Motion
Over the recent years, a variety of reconstructions of the past positions of Australia relative to Antarctica have been put forth (i.e., Cande & Stock, 2004;Whittaker et al., 2007;Williams et al., 2011). We use finite rotations from these studies in order to calculate stage Euler vectors that describe the motion of the Australian plate relative to Antarctica since 70 Ma (Figures 1a and 1b). Toward the beginning of the Cenozoic, the Australian plate was moving west-ward relative to Antarctica at relatively slow rates. This can be inferred from the fact that the stage-rotation pole was located relatively close to Antarctica (see label pre-43 Ma in Figure 1a). Moreover, the stage angular velocity featured values of less than 0.1°/Myr (see Figure 1b). However, the scenario changed after ∼43 Ma when Australia started moving north-ward at faster rates. In fact, its stage angular velocity reached values higher than 0.6°/Myr. Furthermore, the position of the stage-rotation pole moved close to Arabia (see label post-43 Ma in Figure 1a), thus describing a north-ward motion of the Australian plate relative to Antarctica. While these kinematic features are depicted in Figures 1a and 1b, the magnitude of the plate-motion change described above is best evidenced in Figure 1c, where we show the normalized Euler-vector change and the associated 68% confidence range. This illustrates that the largest Australia/Antarctica plate-motion change occurred between 41 and 45 Ma (mid-ages of the two consecutive stage Euler vectors), and is possibly one of the most rapid plate-motion changes ever reconstructed (e.g., Gibbons et al., 2015). Since then, Australia and Antarctica have been spreading apart at relatively fast rates. Figure 1d illustrates, for the period from 70 to 30 Ma, the reconstructed paleo-geography of continents and age of Earth's oceanic lithosphere, which is a proxy for its thickness through a thermal conduction model-for instance through the half-space cooling model (Stein & Stein, 1994). Blue arrows illustrate the motion of the Australian plate relative to Antarctica, upon averaging the available stage Euler vectors over 10 Myr-long intervals. Purple arrows display the absolute motion of the Pacific plate (from the study of Torsvik et al., 2010, which provides Euler vectors over 10 Myr-long intervals). Several tectonics events and observations coincide with the rapid speed-up and change of motion-direction of the Australian plate relative to Antarctica: (i) In the northern hemisphere, the ridge bounding the Pacific plate to the north migrated north-westward between 70 and 50 Ma, eventually plunging into Earth's mantle along the Aleutian trench. This event is also thought to have caused the Pacific plate to change direction of motion, thus generating the Hawaii-Emperor bend at ∼45 Ma (Faccenna et al., 2012;Tarduno et al., 2003;Wessel & Kroenke, 2008;Woodworth & Gordon, 2018). The temporal correlation between these two events, as well as the separation between Australia and Antarctica have been previously noted (Whittaker et al., 2007), but never tested before in the context of models that simultaneously account for torques within the coupled mantle and lithosphere system. (ii) The Indian plate converged towards the larger Eurasian plate at rates among the fastest ever reconstructed (Copley et al., 2010). The final closure of the Tethys Ocean between continental India and China, as well as cessation of spreading between India and Australia caused the merge of these plates into the larger Indo-Australian plate (Liu et al., 1983; see Figure 1d, stages from 50 to 30 Ma). (iii) It has been suggested that temperature anomalies in the lower mantle underneath this region caused the uplift and dynamical sustain of parts of Australia (Czarnota et al., 2013). This, in turn, might have triggered significant pressure-driven flow within the underlying asthenosphere.
We make a preliminary assessment of whether plate boundary forces may provide a plausible explanation for the reconstructed record of Australia/Antarctica separation by estimating their ability to provide the torque-variation needed to generate the reconstructed plate-motion change. Previous studies (Iaffaldano STOTZ ET AL.  Light-color area shows the 68% confidence range (d) Snapshots describing the tectonic evolution of Australia and Antarctica from 70 to 30 Ma, blue arrows show the motion of the Australian plate with respect to Antarctica (this reconstruction), purple arrows are the absolute motion of the Pacific plate (based on Torsvik et al. (2010)). An, Antarctica plate; Au, Australian plate; In, Indian plate; Pa, Pacific plate. Stotz et al., 2018) derived a linear relationship between Euler-vector change experienced by a tectonic plate, and torque-variation necessary to generate a kinematic change. These two vectors are linked by a linear operator that accounts for the shape of the tectonic plate, as well as for the viscosity of the underlying asthenosphere. In Figure 2 we use this relationship and map the distribution of torque-variation magnitudes necessary to change the motions of the Australia plate at 43 Maunder the assumption that Antarctica remains fix during such time period. This distribution has been created from 10 6 samples of the Australian Euler-vector change at 43 Ma-and accounts for an average viscosity of the asthenosphere randomly selected in range from 5 ⋅ 10 19 to 1 ⋅ 10 20 Pa⋅s (e.g., Forte et al., 2010;Lambeck et al., 1996) for each sample of the ensemble. It is evident that the most-recurrent 68% of the ensemble lies in range from ∼1.2 ⋅ 10 26 to ∼2.5 ⋅ 10 26 N⋅m. Based on this argument we can assess under which conditions plate boundary forces may explain the divergence of the Australian plate relative to Antarctica. If one hypothesis that changes of Australian plate boundary forces occurring along a 2,000-4,000 km subducting margins (e.g., Sutherland et al., 2017Sutherland et al., , 2020 are to be responsible for the reconstructed kinematic change, then one would also have conclude that such forces must have changed by 10-20 ⋅ 10 12 N/m, in order to provide the torque-variation needed for the reconstructed plate-motion change. Such a range is indeed too high relative to previous estimates of the typical values of slab-pull or brittle resistive forces along plate margins (e.g., Conrad & Lithgow-Bertelloni, 2002;Faccenna et al., 2012;Stotz et al., 2018). Instead, if one hypothesis that the necessary torque variation arises from changes in the viscous stresses at the base of the Australian plate, then a more plausible average shear stress change in range from 0.4 to 0.6 MPa is needed. The latter range would imply an average pressure-driven Poiseuille flow change of ∼2 cm/yr-i.e., an asthenospheric flow that is faster than the motion of the Australian plate above. Such simple analytical calculations are not intended to provide an explicit test of hypotheses. Nonetheless, they serve well the purpose of providing indications why we focus our numerical-simulations efforts on the above-mentioned hypotheses.

Coupled Models
We use recently developed numerical models, described in detail in Stotz et al. (2017Stotz et al. ( ), (2018, that feature a novel numerical approach and that jointly simulate mantle convection and lithosphere dynamics. Our models build on two previously available codes, one for simulating mantle convection (TERRA), the other for simulating lithosphere dynamics (SHELLS). TERRA is a global, spherical finite element code developed and on which this article is based on are available in Baumgardner (1985) and Bunge and Baumgardner (1995), and further advanced by Bunge et al. (1997); Davies et al. (2013); Stegman et al. (2003); Yang (1997), among others. TERRA solves the classical conservation equations (energy, momentum, and mass) at infinite Prandtl number within a spherical shell domain. Classically, TERRA has been utilized to simulate the socalled mantle circulation mode, in which the distribution of mantle heterogeneity is predominantly controlled by kinematic boundary conditions imposed at the model surface by the user. These typically come from global plate motion reconstructions (e.g., Müller et al., 2008;Seton et al., 2012;Stampfli & Borel, 2002;Torsvik et al., 2010), and are utilized in forward mode (e.g., Bunge et al., 1998;Davies et al., 2012;Müller et al., 2016;Styles et al., 2011) -that is, from sequential data assimilation of past plate motions, rather than adjoint methods (Horbach et al., 2014;Spasojevic et al., 2009). We make use of TERRA in the mantle circulation mode in order to estimate an initial mantle buoyancy field for our joint TERRA-SHELLS simulations at a given starting point in the geologic past-in this study the starting point of our investigations is at 60 Ma. Specifically, we estimate a global mantle buoyancy field STOTZ ET AL.  Histogram for the magnitude of the torque-variation required to generate the Australia/Antarctica plate motion change at 43 Ma. The histogram is created from an ensemble of 10 6 estimates of the torquevariation that account for (i) the uncertainties on the reconstructed Euler vectors and (ii) an average asthenosphere viscosity in range from 5 ⋅ 10 19 to 1 ⋅ 10 20 Pa ⋅ s. that has an inherited 140 Myr-long circulation history, from 200 Ma to 60 Ma, based on the reconstruction of Gibbons et al. (2015). We are aware that other reconstructions have been put forth since the study of Gibbons et al. (2015). However, differences among occur at regional scale, and are negligible at the scale of global plate motions.
We begin by generating global temperature, pressure and velocity fields for the mantle using a modified and benchmarked version of the spherical mantle convection code TERRA. As in the simulations presented in Davies et al. (2012), (2015), calculations are performed on a numerical mesh with ∼80 million discrete nodal points, thus providing the resolution necessary to explore mantle flow at Earth-like convective vigor: Models achieve an internally heated Rayleigh number, based upon reference values, of 5 ⋅ 10 8 , which is similar to estimates of the mantle's Rayleigh number (e.g., Bunge et al., 1997). Our models incorporate compressibility, in the form of the anelastic liquid approximation, with radial reference values represented through a Murnaghan equation of state. Isothermal boundary conditions are prescribed at the surface (300 K) and cosmic microwave background (CMB; 3800 K), with the mantle also heated internally, at roughly chondritic rates (5 ⋅ 10 −12 Wkg −1 ). A free-slip boundary condition is specified at the CMB, whilst surface velocities are assimilated according to plate motion histories from Gibbons et al. (2015), at discrete 1 Myr intervals. Phase changes are incorporated at 410 and 660 km depth, with viscosity varying as a function of depth (z) and temperature (T), following the relation: where T and z are non-dimensionalized by ΔT and mantle depth, respectively, whilst V a and E a are non-dimensional constants controlling the sensitivity of viscosity to depth and temperature. With our choice of V a (2.9957) and E a (4.605), the mean mantle viscosity, which is ∼5 ⋅ 10 20 Pa⋅s in the asthenosphere, increases by a factor of 20 with depth and decreases by a factor of 100 with temperature, providing a reasonable representation of mantle rheology. For initial conditions, a standard convection model is run until a thermal quasi-steady state is achieved. Mesozoic mantle heterogeneity is then approximated by running models with global plate configurations fixed to the oldest available reconstruction at 200 Ma, for ≈100 Myr. Models then evolve from 200 to 60 Myr, when coupling with SHELLS is initiated. Key model parameters are provided in Tables 1 and 2. Thus we begin performing joint simulations using (i) TERRA in order to let the model mantle evolve through geological time and (ii) SHELLS in order to estimates the associated plate motions at the surface. SHELLS is a thin-sheet, finite element code for lithosphere dynamics (see Bird, 1988;Kong & Bird, 1995, for a comprehensive review of the code's technical aspects). Among its features, SHELLS implements the brittle rheology of Earth's lithosphere through input friction coefficients (Iaffaldano, 2012;Kohlstedt et al., 1995;Suppe, 2007). Ductile flow in the lower lithosphere is also taken into account through a parametrized dislocation creep law (Bird, 1998). Lateral variations in the structure of the lithosphere (i.e., Moho depth and thickness of continental roots in continents; depth of lithosphere-asthenosphere thermal transition in oceans) are also incorporated into the finite element description, and constrained from topography and surface heat flow. These are readily available for the present-day, but become more challenging to gather for the geological past. In our simulations we implement the procedure described in Stotz et al. (2017) for incorporating first-order past lithospheric structures into our model finite element grids, which are made available as part of this study. From these geological constrains, as well as from the asthenospheric flow STOTZ ET AL.  calculated in TERRA, SHELLS is able to estimate all torques acting upon tectonic plates at a given point in the model past (i.e., at 60 Ma in the cases explored here), and to output global plate motions. Next, these are passed onto TERRA and used as a surface kinematic boundary condition in order to evolve the model mantle flow through geological time. In this process, we also utilize the lateral thickness variations cast into the finite element lithosphere grid of SHELLS in order to modify the lateral viscosity variations cast into TER-RA's finite element grid. This allows us to include a viscosity structure in the upper layers of TERRA that is consistent with the one of SHELLS. By iterating this procedure through model time, our joint TERRA/ SHELLS models output absolute global plate motions that are then converted into Euler vectors for each plate, in a post-processing step. This allows us to calculate model Euler vectors at a temporal resolution of 1 Myr that can be compared against the available reconstructions.

Case Studies
We examine the tectonic evolution of the Australia/Antarctica region in four different scenarios, in order to test hypotheses on the mechanisms that might have led to the separation between the Australian and Antarctica plates at ∼45 Ma. To this end, we follow a similar methodology presented in Stotz et al. (2017): we build four sets of SHELLS finite element grids that cast the global lithosphere structure at a temporal resolution of 1 Myr. These four sets, however, describe different evolution of the lithospheric and mantle buoyancy fields. They are depicted in Figure 3 and described in the following: (i) The first set of grids follows the global tectonic evolution put forth by Gibbons et al. (2015) (see Figures 3a-3c). Specifically, the Pacific ridge subducts underneath the Aleutians at ∼50 Ma, and the India/Australia ridge (known at present-day as the 90 East Ridge) heals after ∼45 Ma. Furthermore, the mantle buoyancy field (shown in Figures 3d and 3e) is drawn from the mantle circulation history modeled after imposing a kinematic surface boundary condition from 200 Ma to 60 Ma-i.e., it features a 140 Myr-long inherited mantle circulation history. Instead, from 60 Ma towards present-day global surface velocities are estimated by SHELLS through the procedure described above. We label this as our Reference Simulation (REFS). The mantle temperature distribution at 60 Ma for 100 and 2,300 km depth is shown in Figures 3d and 3e. Note that at 2,300 km depth one can recognize warmer-than-average areas of upwellings that are related to the African and Pacific super swells-these are also referred to in the literature as large low-shear-velocity provinces (LLSVPs). (ii) In order to test the effect that the subducted Pacific ridge might have had on the separation of the Australian and Antarctica plates, we build a different set of grids that do not include such a subduction event. Figures 3f-3h displays finite element grids where the Pacific ridge remains unsubducted after 50 Ma. In other words, the Pacific plate is deliberately not allowed to subduct along the Aleutians-Japan trenches. Since in these models, the ridge never subducts, one would expect that the Pacific plate will not be able to change direction of motion. We refer to this simulation as the No Subduction of the Pacific Ridge (NSPR). In this particular simulation, the mantle buoyancy is the same of REFS (see Figures 3i and 3j)  (iv) Lastly, we perform a simulation that features a random distribution of mantle buoyancy at the initial model time of 60 Ma. This is obtained using a random mantle temperature distribution as initial condition, then running a mantle circulation simulation for 50 Myr of model time with the kinematic boundary condition at the surface kept systematically fixed to the reconstructed plate motions at 60 Ma. This is a computationally necessary step to take in order to generate a shallow, sub-lithospheric buoyancy field that, while not being representative of the actual history of Earth's mantle circulation, may still be used to calculate slab-pull forces at subducting margins globally. We refer to this simulations as the No Inherited Mantle History (NIMH). Its purpose is to obtain a mantle buoyancy field where sinking lithospheric slabs only reach upper-mantle depths of ∼1,300 km, without however penetrating into the lower mantle. Therefore, in this simulation the Poiseuille component of the asthenospheric flow is not representative of the past history of mantle circulation-although it still complies with the laws of physics cast into the conservation equations mentioned above. Figure 3s shows the temperature field at 100 km depth, where slabs are visible along convergent margins. As expected, the temperature field is, to a certain degree, similar to the one REFS (see Figures 3d and Figure 3s). However, STOTZ ET AL.

Results
We use the output Euler vectors for the absolute motions of the Australian and Antarctica plates to calculate Euler vectors for their relative motion at a model temporal resolution of 1 Myr. Next, we use these model Euler vectors to estimate the Euler-vector change through time. We do so by first interpolating the model Euler vectors at the same temporal resolution of the relative motion reconstructed from the data of Cande and Stock (2004); Whittaker et al. (2007); Williams et al. (2011), and then taking the Euler-vector difference between each time period. This allows us to mimic the temporal resolution of the Euler vectors reconstructed from actual pickings of the ocean-floor magnetization pattern (e.g., Seton et al., 2014). In Figure 4 we assess which of our numerical simulations yields kinematics that are comparable with the reconstructed Euler-vector change that describes the separation of Australia and Antarctica. In blue, we show the normalized Euler-vector change between ∼60 and ∼20 Ma arising from kinematic reconstructions, along with the associated 68% confidence range. Simulations REFS (in dashed black in Figure 4) reproduces well the reconstructed Euler-vector change. On the one hand, this means that such a simulation allows capturing the timing and magnitude of the separation between the Australian and Antarctica plates. On the other hand, however, such a result alone does not allow us to single out the specific process responsible for the Australia/Antarctica separation, since in REFS several processes (described above) are acting simultaneously.
Therefore, we evaluate the impact that subduction of the Pacific ridge prior to ∼50 Ma has on the separation of Australia and Antarctica by calculating the model Euler-vector change for simulation NSPR (in magenta in Figure 4). The results of this simulation demonstrate that the subduction of the Pacific ridge played virtually no role in the Australia/Antarctica separation, since removing such a process from the modeled dynamics yields a temporal pattern of Euler-vector change very similar to the reconstructed one. In other words, there is no modeling basis to argue that the temporal correlation between these two tectonic events is in fact evidence of a causal relationship. We interpret this in light of the fact that the margin between the Pacific and Australian plates comprises several small plates that may act as a buffer and prevent lithospheric stresses from being fully transmitted from the Pacific realm onto the Australian plate. Our simulations further indicate that the merging of the Australian and Indian plates into the Indo-Australian unit did not play a role in the separation of Australia from Antarctica (i.e., NHIR). This is evident from the predicted Euler-vector change pattern (dashed green in Figure 4). As in the previous test, the Euler-vector change pattern displays a maximum variation at the time of Australia/Antarctica separation. In other words, our results suggest that the merge between the Indian and Australian plates is not sufficient to generate enough pulling force to drag along the Australian plate northward. We concede that our models do not allow for the possibility to incorporate any diffuse deformation region separating plates, as, for instance, the one developed within the Indian Ocean region (e.g., Gordon, 1998;Gordon et al., 1998;Krishna et al., 2009;Wiens et al., 1985). However, we note that including such a deformation within the newly formed Indo-Australian unit would only hinder the ability to pull the Australian plate away from Antarctica.
Instead, the Euler-vector change pattern arising from NIMH is significantly different from the reconstructed one. In the NIMH simulation, the initial mantle buoyancy field at 60 Ma does not represent the mantle circulation history between 200 and 60 Ma. Arguably, the buoyancy field in the lower mantle and the locations of subducting slabs control to a large degree the asthenospheric Poiseuille flow (e.g., Cande & Stegman, 2011;Morgan et al., 1995;Stotz et al., 2018). In fact, hot buoyant material rises upwards from the lower mantle and reaches the asthenosphere, thus triggering positive pressure-driven flow. At the same time, cold downwellings also contribute to the pressure-driven flow within the asthenosphere in the form of negative pressure signals. Because of this, the Poiseuille component of the asthenospheric flow in the NIMH simulation, while disjointed from the actual circulation history of Earth's mantle between 200 and 60 Ma, remains representative of the upper-mantle buoyancy field associated with the 60 Ma lithosphere kinematics. For this reason, we performed two distinct simulations of joint mantle and lithosphere dynamics associated with the NIMH buoyancy field. In one, lithospheric plates are subject to both the slab-pull forces associated with the downwellings of the NIMH buoyancy field and the shear stresses associated with the NIMH asthenospheric flow. The normalized Euler-vector changes from this simulation is the upper dashed red line in Figure 4. In the other simulation, lithospheric plates are subject only to the slab-pull forces associated with the downwellings of the NIMH buoyancy field. We achieved this by imposing a passively flowing asthenosphere (i.e., one where only the Couette component of flow is kept) in the NIMH convection model, which has been derived using a random mantle temperature distribution as initial condition (see details above). The normalized Euler-vector changes from this simulation is the lower dashed red line in Figure 4. These two simulations represents the end-members of a range of scenarios where the history of mantle circulation has virtually little to no impact on the dynamics of plates. The average between the two end-member simulations is the solid red line in Figure 4. When the mantle buoyancy history is not inherited, the normalized change of the Australia/Antarctica Euler-vector decreases by 30%-80% (light-red area in Figure 4), falling outside the 68% confidence range warranted by kinematics reconstructions.

Discussion
Taken together, our simulation results support the notion that the active asthenospheric flow underneath Australia and Antarctica is ultimately responsible for their separation, since removing it from our numerical models yields a poor fit to kinematic reconstructions. We find useful to discuss such flow in relation to two main, well-established characteristics of global viscous convection patterns of the Earth's mantle. The first, the mantle buoyancy field at a given time is shaped to a significant degree by the history of subducting slabs, which is illustrated in global collections of past plate-motion reconstructions (e.g., Gibbons et al., 2015;Müller et al., 2008;Seton et al., 2012;Stampfli & Borel, 2002;Torsvik et al., 2010). When used as a surface boundary condition in mantle circulations models, these reconstructions result in lateral variations of mantle temperature that are in line with seismic tomography models (Grand et al., 1997;Davies et al., 2012Davies et al., , 2015Ritsema et al., 2011;Romanowicz, 2003;Schuberth et al., 2009). Common to these models is the fact that the lower mantle exhibits two very distinct temperature anomalies at the core mantle boundary: One below the Pacific basin, while other underneath the African continent. Furthermore, these features are bounded by the cold down-going slabs in correspondence at subduction zones -e.g., the subduction zones along the entire North and South American continents. This defines a spherical degree-two pattern of mantle convection cells. In simple words, it means that Earth's mantle features two main areas of upwelling encapsulated by cold down going slabs. Indeed such a characteristic feature can be seen in our REFS (Reference simulation), which account for the 140 Myr-long history of global subduction. Specifically, a Figure 5a displays global temperature anomalies in the lower mantle (2,300 km depth) at 60 Ma of model time, when coupling of mantle and lithosphere dynamics is initiated. In contrast, the pattern of lower mantle temperature anomalies does not feature an equally strong degree-two component when inheritance of the 140 Myr-long history of global subduction is suppressed, as in the NIMH simulation (see Figure 5b).
The recently proposed Cathless parameter (Hoeink & Lenardic, 2010;Richards & Lenardic, 2018) is helpful in understanding how the long term evolution of the mantle flow influences plate dynamics. The Cathless parameter is puts in relationship the viscosity contrast between upper mantle and asthenophere, as well as the ratio between the length of a convection cell and the depth of the mantle. Its power lies in enabling a prediction of whether a pattern of mantle convection falls in the Poiseuille-(mostly pressure-driven) or in the Couette-driven (mostly driven by upper plate motions) regime. For a degree-two convection (aspect ratio of ∼4) and a viscosity contrast of 10 2 the Cathless parameter implies a Poiseuille-driven flow (see Figure  10 in Richards & Lenardic, 2018). This indicates that in the REFS simulation (associated with the pattern in Figure 5a) plate motions at the surface are largely influenced by the pressure-driven Poiseuille flow within the asthenosphere. This is in line with studies on the dynamics of the Pacific plate (e.g., Stotz et al., 2018). In contrast, this is not the case for the NIMH simulation.
The second characteristic is that previous studies (Bunge et al., 1998) showed the typical time-scales over which the large-scale (i.e., scale of large tectonic plates) pattern of mantle convection varies is on the order of 100 Myrs or more. This means that a given pattern of global convection, once established, it keeps its first-order-degree features for such long periods of time. In the context of the hypotheses tested here, this means that the pattern of convection resulting from the pre-60 Ma circulation mantle history and responsible for the Australia-Antarctica separation around 43 Ma will still be driving their present-day divergent motion. In line with the arguments above, the REFS simulation (which features a significant degree-two component of flow) exhibits fast flow velocities within the asthenospheric channel. In Figure 6a we show the model mid-asthenospheric temperature and flow fields in the region comprising the Indian ocean, the Australian and Antarctica plates, and the Pacific basin for the REFS simulation. In particular, in our models the lithosphere thickness between Australia and Antarctica is constrained by the ocean floor age. Thus, the rift zone in our models is compromised by thin ocean floor and bounded by two thicker continental regions. In our models, the hot temperature anomalies are not necessarily located directly underneath the Australia and Antarctica rift margin, in fact they come from the east and then flow north underneath the Australian continent. Instead, Figure 6b illustrates flow in the same region, but for the NIMH simulation, where the inherited mantle circulation history from 200 to 60 Ma has been removed so that only the downwelling buoyancy field associated with subducted slabs generates lateral pressure gradients.
A comparison of the two panels illustrates how the inherited mantle circulation history controls asthenospheric flow, and how it drives the motions of tectonic plates through shear stresses at the lithosphere base. Specifically, Figure 6a displays much stronger temperature anomalies within the asthenosphere relative to Figure 6b. For instance, two upwellings featuring high temperature anomalies of about 400 K are clearly visible in the Indian and Pacific basins. In contrast, Figure 6b displays much weaker upwellings featuring temperature anomalies of 200 K or less. Furthermore, the pattern of asthenospheric flow in Figure 6a is stronger and more coherent over longer distances than it is in Figure 6b, as for the northward flow underneath Australia or the westward flow underneath the Pacific basin in Figure 6a. Flow velocities there are ∼10 cm/yr or higher, while in Figure 6b flow velocities are systematically less than ∼7 cm/yr. These strong differences arise because in the NIMH simulation, the pattern of temperature in the lower mantle originates from a random initial mantle structure. Consequently, while the Couette-type flow imposed onto the asthenosphere by subduction generates regions of negative pressure and thus fast flow near subducting slabs (e.g., Aleutians-Japan convergent margin), the randomness of the initial mantle structure overall results in weak positive pressure regions within the asthenosphere. This in turn causes a weaker flow field as well as temperature variations occurring over relatively short distances, and is responsible for generally smaller convection cells (see Figure 3t) that characterize the NIMH simulation. Furthermore, low temperature anomalies are not capable of generating large-enough topographic swells (i.e., dynamic topography), as instead high temperature anomalies are. In other words, the distribution of the length scale temperature anomalies in our reference simulation provides enough buoyancy over greater distances and, thus, is capable of generating dynamic topography signals within the Australian plate, such characteristics are in agreement with the patterns Cenozoic dynamic topography recorded in and around the Australian continent (Czarnota et al., 2013(Czarnota et al., , 2014.

Conclusions
In this study, we analyzed the mechanisms that may have lead to the separation of the Australian and Antarctica plates after ∼45 Ma. To do so, we collected finite rotations available in the literature to estimated stage Euler vectors that describe the relative motion of the Australian plate with respect to Antarctica (Can- de & Stock, 2004;Whittaker et al., 2007;Williams et al., 2011). Next, we estimated the temporal Euler-vector changes to determine the timing of the separation between Australia and Antarctica. As several tectonic events coincide with the timing of maximum Euler-vector change, we tested through global models of the coupled mantle and lithosphere dynamics which of these events is capable of driving the separation of the Australian and Antarctica plates. Specifically, we analyzed three plausible mechanisms that might have led to the separation of the Australian plate from Antarctica. We found that the subduction of the Pacific ridge after ∼50 Ma and the consequent change in direction of the Pacific plate at the time of the Hawaii-Emperor bend played virtually no role in separating Australia and Antarctica. Similarly, the healing of the India/Australia at ∼43 Ma, merging Indian and Australian plates into the Indo-Australian unit, had no clear impact on the separation of the Australian plate from Antarctica. Our results, however, indicate that the separation between the Australian and Antarctica plate was driven to a large extent by pressure-driven asthenospheric Poiseuille flow associated with the mantle buoyancy field that is inherited from the Mesozoic circulation history.