Constraining the Range and Variation of Lithospheric Net Rotation Using Geodynamic Modeling

Abstract Lithospheric net rotation (LNR) is the movement of the lithosphere as a solid body with respect to the mantle. Separating the signal of LNR from plate tectonic motion is therefore an important factor in producing absolute plate motion models. Net rotation is difficult to constrain because of uncertainties in geological data and outstanding questions about the stability of the mantle plumes used as a reference frame. We use mantle convection simulations to investigate the controlling factors for the magnitude of LNR and to find the statistical predictability of LNR in a fully self‐consistent convective system. We find that high lateral viscosity variations are required to produce Earth‐like values of LNR. When the temperature dependence of viscosity is lower, and therefore slabs are softer, other factors such as the presence of continents and a viscosity gradient at the transition zone are also important for determining the magnitude of net rotation. We find that, as an emergent property of the chaotic mantle convection system, the evolution of LNR is too complicated to predict in our models. However, we find that the range of LNR within the simulations follows a Gaussian distribution, with a correlation time of 5 Myr. The LNR from the models needs to be sampled for around 50 Myr to produce a fully Gaussian distribution. This implies, that within the time frames considered for absolute plate motion reconstructions, LNR can be treated as a Gaussian variable. This provides a new geodynamic constraint for absolute plate motion reconstructions.

Absolute plate motion reconstruction models all produce some LNR, although the magnitude varies depending on the reconstruction. This variation is partly due to differences in the choice of reference frame used. The creation of absolute plate motion models requires high-quality, widespread palaeomagnetic data which is only available for the relatively recent geological past. The further back in geological time one goes, the sparser and more unreliable the palaeomagnetic data becomes, reducing the certainty with which absolute plate motions can be reconstructed, and therefore estimates of LNR also become more unreliable.
The range of LNR values predicted by absolute plate motion models goes from around 0.13°/Myr (Torsvik et al., 2010) to 0.44°/Myr (Gripp & Gordon, 2002) when a deep hotspot reference frame is used, potentially increasing to  E 1°/Myr if shallow reference frames are used (Doglioni et al., 2015), or as low as 0.11°/Myr when the models are setup to explicitly minimize LNR (Tetley et al., 2019). Studies of seismic anisotropy, which can be used as a proxy for shearing between the lithosphere and the mantle suggest amplitudes of 0.2-0.3°/Myr (Becker et al., 2008(Becker et al., , 2014Conrad & Behn, 2010;Kreemer, 2009), which are within the range proposed by plate models. In the absence of constraints on LNR, one solution is to minimize LNR at all times in plate motion reconstructions (Tetley et al., 2019). However, these models generally aim to fit geological and palaeomagnetic data but are not constrained by the physical laws that control geodynamic flow.
In order to improve absolute plate reconstructions and extend them further into the past, we need to know how much and over what time scales LNR can reasonably vary. Geodynamic modeling lets us explore the ranges and timescales of variation, allowing us to assess whether the LNR predicted by plate reconstruction models is geophysically reasonable. Ideally, we would also be able to identify tectonic situations which lead to high or low net rotation. We could then make an estimate of the likely LNR in a particular global tectonic configuration even if there is insufficient palaeomagnetic data to constrain a reconstruction well. In order to do this, we investigate the origins of LNR in mantle convection models.
The movements of the convecting mantle can be decomposed into poloidal motions, the up-and downwellings caused by buoyancy differences, and toroidal motions, the origin of which is generally attributed to lateral viscosity variation (Bercovici et al., 2000). Net rotation of the lithosphere is the degree 1 of the toroidal component of velocity. In this study, we use two-dimensional simulations. The poloidal motions and degree 1 toroidal component (LNR) can be studied in 2D, depending on the geometry of the model. Although we lose the higher orders of the toroidal component by using 2D models, the computational saving compared to using expensive 3D models means we can investigate a greater range of parameter combinations to constrain some basic relationships between convection parameters and net rotation. We use the mantle convection code StagYY (Tackley, 2008) with spherical annulus geometry (Hernlund & Tackley, 2008). Unlike in 2D models which use an axisymmetric geometry, net rotation of the lithosphere is possible in a spherical annulus. Similar studies have also been performed using a 2D cylindrical geometry (Gérault et al., 2012), but spherical annulus geometry has the advantage that velocities are closer to those expected in a full 3D simulation (Hernlund & Tackley, 2008).
The convection models allow us to investigate which parameters control the magnitude and variability of LNR. It has long been recognized that lateral viscosity variations are required to produce net rotation (Alisic et al., 2012;Becker, 2006;Gérault et al., 2012;Ricard et al., 1991;Zhong, 2001), Continents provide a strong lateral viscosity variation (Becker, 2006;Rudolph & Zhong, 2014;Zhong, 2001), enhanced if they have deep roots, as do temperature and stress-dependent viscosity laws (Alisic et al., 2012). Deep mantle viscosity structure affects the speed of plate motion, by determining the sinking rate of slabs. A viscosity increase at the 660 km transition zone may slow or stop slab penetration into the deep mantle (Conrad & Lithgow-Bertelloni, 2002), while increasing temperature dependence of viscosity in the lower mantle increases the viscous drag around the slab tip, slowing its sinking (Alisic et al., 2012). Other factors that have been proposed to affect net rotation include the position of oceanic ridges, trench migration, slab properties, and a weak asthenospheric layer. Gérault et al. (2012) investigated these factors in 2D convection models and found that asymmetric slab dips and ridge positions promoted net rotation. However, geodynamical studies tend to find NR rates which lie at the lower end of the ranges proposed by plate tectonic models, generally less than 0.15°/Myr (Becker, 2006;Ricard et al., 1991;Zhong, 2001), except the models of Alisic et al. (2012), which have NR values up to around 0.25°/Myr. The earlier geodynamic modeling studies did not use highly temperature-dependent viscosity. Slabs in the mantle could therefore be relatively soft. In models that have larger lateral viscosity variations, such as those used in this study and Alisic et al. (2012), slabs remain stiff, and values achieved for net rotation are in line with those inferred for Earth.
In this work, we study net rotation in evolving mantle models. Unlike previous studies, we have both high lateral viscosity variation, and therefore stiff slabs in the mantle, and continuously evolving mantle structures. Our models are neither initialized nor constrained by the imposition of present day observations of the surface or mantle structure. This allows us to investigate how net rotation evolves as a statistical and emergent property of a fully self-consistent mantle system, where all tectonic behavior arises naturally from physical laws, rather than relying on the imposition of surface forcing. We can therefore study the constraints on the magnitude and statistical predictability of LNR in our models, independent of uncertainties and inconsistencies arising from the geological record.
We perform a parameter space search in order to more fully investigate the effects of different geodynamical properties on the magnitude and variation of net rotation. This gives us a statistical way in which we can constrain the net rotation of the Earth's lithosphere for use in absolute plate motion models. We vary a variety of convection parameters that have been identified in previous studies to have effects on the magnitude of net rotation. The models we use can include continents (although not all of our models do) and they develop self-organizing plate tectonics with spreading ridges and subduction zones. In agreement with previous studies, we find that changing the viscosity, especially the parameters that increase the lateral variation of viscosity, increases the magnitude of lithospheric net rotation. For models with highly temperature-dependent viscosity and therefore stiff slabs, we achieve ranges of LNR within the published ranges for Earth. Our large viscosity contrasts mean that changing the presence and thickness of continents makes little difference to LNR, because the presence of stiff slabs in the mantle creates a much greater viscosity contrast than the continents, in contrast to previous studies.

Methods
Our convection models are incompressible and were run using the geodynamical code StagYY with a 2D spherical annulus geometry (Hernlund & Tackley, 2008). The resolution is 1,152 cells azimuthally by 128 radially, with  6 1.6 10 E Lagrangian tracers. The surface and core-mantle boundaries both have free-slip conditions. Where present, two continents are initialized to cover 30% of the surface. They are stiffer and more buoyant than the upper mantle and lithosphere to ensure that they remain at the surface. They have lower yield stress margins, allowing them to form a single supercontinent and subsequently break up. The mantle is initialized with a layer of denser primordial material at the core-mantle boundary. There is internal heating, but this does not decay with time. We, therefore run the models until they reach a stable whole mantle mean temperature, at which point we judge them to have entered a statistically stable convective state. Once the models have reached a statistical steady state, we run them for at least a further 600 time steps. We use this stable phase of the run for all of our subsequent analyses.
We use a Rayleigh number (   there is no variation in this axis. The velocity and gravity vectors also have no components in this neglected dimension. For net rotation, this means that only degree 1 toroidal motion is possible because higher degrees require variations in velocity perpendicular to the slice. Unlike when modeling the mantle using an axisymmetrical geometry, degree 1 net rotation is possible in spherical annulus geometry because it is a laterally uniform equatorial slice. In contrast, the axisymmetric system has internal impermeable reflecting boundaries at the poles, preventing any toroidal velocity component.
In order to solve the Navier-Stokes equations, it is necessary to define a reference frame for velocity. There are different ways to define the reference frame. StagYY uses a no surface net rotation reference frame to calculate the velocity at each point. Technically speaking, a numerical solution is obtained in an arbitrary reference frame, and a global solid rotation is imposed to move this solution back to the no surface net rotation reference frame. We then post-process the models to find surface velocity and rotation in a no mean

Viscosity Activation Volume
An increase in viscosity activation volume, a E V (Equation 2), generally, but not always, leads to an increase in LNR. The increase can be up to 1 order of magnitude for a 5-fold increase in a E V , although for most parameter combinations, a 3-5-fold increase in a E V produces around a doubling of LNR. The only sibling set where a E V does not increase LNR has very thick continents and a E E = 1. The two models have LNR with the same order of magnitude, therefore we attribute the decrease in LNR with increased a E V to the chaotic nature of the evolution of LNR and suggest that the models are nearly statistically identical. Increasing a E V also increases temperature dependence of viscosity, although less directly than a E E (see Equation 2). It has a correspondingly smaller effect on LNR, with magnitude variation independent of a E E .

Viscosity Jump
A factor 30 increase in viscosity at the transition zone almost always produces an increase in LNR. The implementation of the viscosity jump in our models effectively reduces the viscosity of the upper mantle compared to models without a jump, whilst the lower mantle has comparable viscosity in all models. The addition of a viscosity jump increases LNR between comparable models by one to two orders of magnitude, with a greater increase in models with lower a E E . In models with higher a E E , a reduction in viscosity by a factor 30 produces a much less pronounced increase in lateral viscosity variations because the thermal effects on viscosity are much greater than  E 30.

Continent Presence and Thickness
At low values of a E E , adding continents increases LNR by up to an order of magnitude. However, with increasing a E E , this relationship breaks down. The thickness of the continents, once they are present (thickness  E 0), makes very little difference to LNR at any a E E . Previous studies found that both the presence and thickness of continents (Becker, 2006;Zhong, 2001) had an impact on LNR rates. However, these studies have much lower lateral viscosity variations (e.g., lateral viscosity variation of around two orders of magnitude in Becker, 2006), whilst ours have lateral viscosity variations of over seven orders of magnitude. When slabs  can produce such a large lateral viscosity variation in the mid-mantle, the smaller viscosity variation produced by stiff continents at the surface is of secondary importance. Our low a E E models show similar results to previous studies, because they share similar rheological properties.
No other parameter makes a reliable difference to LNR. Very high yield stress causes the models to stagnate. This tends to lead models to have relative time-invariant values of net rotation. In some cases, the LNR remains around zero, but it can also lead to moderately high magnitude uni-directional net rotation. These simulations, therefore, appear to have higher magnitudes of LNR/ V rms E when the stable runtime average is plotted in Figure 1, but are in fact relatively inactive compared to the other models. Changing the thickness and viscosity of the layer of dense material at the base of the mantle makes no difference to the LNR. Increasing the internal heating through the decay of radioactive elements also makes no difference to the magnitude LNR, despite this affecting the convective vigor. Variation within each sibling set for these parameters is due to the chaotic nature of mantle convection, meaning that a small variation can change the convection patterns. With longer run times, we expect these model sets to converge to have stable, near-identical average rates of LNR.

Viscosity and LNR
In Section 3.1, four convection parameters were identified that had a discernible effect on lithospheric net rotation. Three of these parameters directly control the viscosity structure of the models. The fourth, continental lithosphere, introduces large viscosity variations at the surface. The relationship between viscosity and LNR becomes much more obvious when just the upper mantle is considered ( Figure 2). In the lithospheric upper mantle (at depth 0.04, equivalent to 115 km), we see a log-linear relationship between the minimum viscosity at each time step, averaged over the models' run, and the root mean squared LNR. We use the minimum viscosity because at 115 km the continental roots all have the same maximum viscosity. The continents, therefore bias the mean viscosity, meaning that models with and without continents cannot be compared. We show the same results twice, first with LNR and then with the LNR scaled by the RMS of the surface velocity.
The presence of continents makes no significant difference to the relationship between viscosity and LNR. The relationship is dependent on the minimum viscosity under the oceans, with the relationship being almost identical when the continental lithosphere is excluded (middle and right columns, Figure 2) because the minimum viscosity almost always occurs under the oceans. The minimum viscosity, and therefore the net rotation, is dependent on the viscosity parameters chosen in the model, namely the presence of a viscosity jump and the activation energy of viscosity, which determines the sensitivity of viscosity to temperature variations. However, when LNR is scaled by Vrms, the relationship weakens. This is expected because much of LNR magnitude is determined by the vigor of convection, which is a function of viscosity.

LNR Variation With Time
The results presented so far were all for the RMS value of net rotation over the runtime of the model, with some correlations between parameters and LNR found when we consider the full suite of models. Whilst this gives a statistical average, it does not give details about the expected value at any time step nor how the rotation varies within the model. For absolute plate motion models, estimates of how much LNR varies with time are required. We therefore investigate if particular tectonic and convective configurations can be linked to LNR magnitude and the timescales over which LNR varies.
Visual inspection of the evolution of individual models does not yield any clear relationship between geodynamic patterns and changes in net rotation. Several models change from having regular supercontinent accretion and break-up cycles to having permanent supercontinents. In most of these models, but not all, the phase where the continents regularly disperse and converge is associated with lower LNR. The RMS surface velocity does not change when the models switch to a single, stable supercontinent, therefore the increase LNR is not caused by an increase in convective vigor when a single supercontinent dominates. It seems that the frequent change in direction of continent motion as they converge and dispersal constraints the LNR, probably by maintaining higher degrees of symmetry in the model. However, during this alternating phase, the net rotation does not increase when the continents have come together.
The majority of models do not have supercontinent cycles. The continents converge during the pre-stabilization phase before analysis starts and do not subsequently diverge. In these models, there is no clear link between geodynamic processes and LNR. We present model 12 (see Appendix A for details) as a representative model, having mid-range viscosity parameters (   3.0 a a E E V ) and continents. The maximum LNR is 0.184°/Myr, although for most of its life the LNR does not exceed 0.05°/Myr. Figure 3 shows the variation of several convection metrics. We see that variations in temperature and convective vigor are not clearly correlated with LNR. LNR varies on much shorter timescales than temperature. Surface velocity also varies on shorter timescales than temperature, but peaks in surface velocity are not correlated with peaks in LNR.
Sixteen snapshots of the temperature of the mantle from model 12 are shown in Figure 4. They are ordered according to LNR, with low, medium, and high values, and samples from the single very large peak in LNR at 4,000 Myr. The snapshots in the first three groups are all far enough apart in age to be uncorrelated. There is no clear difference between snapshots with low or high LNR. They all have active subduction zones with slabs extending to the lower mantle, with at least one in each category having a slab adjacent to a continent. All have at least three ridges with the hot, thin oceanic lithosphere, and they all have upwellings and slab remnants in the lower mantle. The last row are snapshots taken from the peak of net rotation at 4,000 Myr, which is significantly above average for the model. At the beginning of the peak (3,959 Myr), both slabs have broken off, disconnecting the lithosphere from the mantle. This could explain why the lithosphere is free to rotate. However, we see this elsewhere in the model's history without a corresponding peak in LNR. The cold slab remnants in the lower mantle are also asymmetrically distributed, but we see similar uneven distribution in the top left snapshot with very low LNR. What has caused this peak is therefore unclear. This model is representative in demonstrating the difficulty of finding a relationship between mantle structure and LNR. We have investigated a variety of metrics from the models, using viscosity, temperature, composition, radial, and azimuthal velocity, to which we have applied various image and signal processing methods, but we have been unable to find a robust and reliable predictor of temporal variations in LNR.
Because we are unable to find a clear relationship between mantle structure and LNR by simply inspecting our models, we turn to statistical descriptions of mantle structure in order to describe how LNR changes as models convect. This allows us to make generalized descriptions of the models' behavior, even if we cannot say precisely what the LNR of the models will be at a given time. is dependent on the minimum asthenospheric viscosity and the surface velocity from Figure 2, we can estimate the range over which LNR varies within our models.
We then investigate how the 1D profiles of temperature and viscosity are correlated with LNR during the stable runtime of each individual model. We can see if particular changes within a model tend to be linked to a change in net rotation, for example, a decrease in temperature, generally caused by increased subduction. This gives us a method to identify convection features that may control LNR, even if they are not clear from inspection of the models. We use the Spearman rank correlation coefficient to correlate changes in the 1D mantle structure during the stable runtime of our models with changes in LNR. This measures the monotonic relationship between two parameters and has the advantage that the relationship does not have to be linear in order to produce a high correlation. The Spearman rank correlation coefficient is simply the Pearson correlation coefficient between the rank variables of two parameters, E x and E y , and can be calculated by: where i E rx and i E ry are the ranks of E x and E y values respectively, ranging from 1 to E n .
The correlation between the magnitude of LNR and 1D profiles of temperature and viscosity is calculated separately for each individual model, and the stacked results for all the models are shown in Figure 5. Positive correlations mean that magnitude of LNR increases with the parameter in question, negative that one increases as the other decreases. There is no feature that is universally correlated with LNR across all models. However, when individual models are considered, some relevant correlations emerge.

Correlation of Thermal Profile With LNR
Each set of 1D correlation profiles are stacked according to the convection properties of the models, with increasing values of a E E and a E V in the first two columns and then the models are separated depending on if they have a viscosity jump in the last two columns. The first row shows the correlation between thermal structure and LNR for each model. The correlation between thermal structure and LNR varies only slightly depending on the convection parameters used. Near the surface, the thermal structure is positively correlated with LNR for depths above 0.096, which corresponds to the base of the continents in the majority of models. The correlation is weak, below 0.5 in all cases. Going deeper, the upper mantle temperature, below the lithosphere, is controlled by the changing balance between cold down-going slabs and hot upwelling material. Colder upper mantle has a weak correlation with increased LNR because a more mobile surface produces more slabs which cool the upper mantle. Below this, no pattern in correlation emerges. Models with low a E E and a E V show very little correlation between thermal structure and LNR, although this is not surprising since these simulations have very low LNR in the first place. Simulations without a viscosity jump show a greater change in correlation direction below the continents because they have a higher upper mantle viscosity, making them more temperature-sensitive. In simulations with a transition zone viscosity jump, fewer than half show any change in correlation around the transition zone (depth 0.286).

Correlation of Viscosity Profile With LNR
The strength of the correlation between the 1D viscosity profile and LNR is very clearly a function of the viscosity parameters chosen. Increasing a E E and a E V both increase the positive correlation in the mid-mantle.
These parameters enhance the viscosity contrast of slabs, highlighting that at times when there are more slabs in the mid-mantle, net rotation is often stronger. The effects of a fast-moving asthenosphere are not visible in the viscosity correlation. This is because the sensitivity of viscosity to temperature is somewhat lower at lower pressures. The correlation is also stronger in simulations without a viscosity jump, because these have higher upper mantle viscosity. We also consider the viscosity profile below the oceans and continents separately. Where present, the continents cover 30% of the surface. The sub-oceanic mantle makes up the majority and it is therefore unsurprising that the sub-oceanic correlations are more similar to the global viscosity correlations. The sub-continental mantle shows very little correlation with NR in any simulation, although we do see a slight increase in correlation immediately below the continents (depth 0.096) in some simulations. This is again related to surface mobility -stationary continents insulate the upper mantle, reducing viscosity. More mobile continents move over colder areas and so will have slightly higher subcontinental viscosity, which correlates with higher mobility and net rotation. The positive mid-mantle correlations in high a E E and a E V cases continues for sub-oceanic profiles. Most slabs are subducted in sub-oceanic mantle so this is to be expected. The effects of the continents can be seen with an increase in correlation strength at a depth of 0.096. Above this depth, the contribution of continents limits the correlation between viscosity and LNR.
Whilst there is no clear pattern linking thermal or tectonic structure to LNR, we can at least investigate over what time scales LNR is stable. Figure 6 shows a cumulative histogram of the change in LNR after varying time lapses. When the difference between two random Gaussian signals with a standard deviation of 1 is calculated, 60% of samples would be within approximately 0.75 of each other. Figure 6 suggests that after a time lapse of 100 Myr, the correlation between two measures of LNR is no better than random. With smaller  We also investigate over what timescales we need to sample net rotation in order to produce a Gaussian distribution. This is an important consideration because it determines over what timescales we can reasonably use a Gaussian distribution to model LNR. In order to do this, we separate the LNR for the whole model into chunks of set length. We repeat the exercise, with chunk lengths between 10 samples up to 1,000 samples. We then normalize the LNR in each chunk. If the distribution within the chunk is Gaussian, the normalized chunk distribution will be Gaussian. To test this, we use the Kullback-Leibler ( KL E D ) divergence between the chunk distribution and a Gaussian distribution. If the distribution of LNR with the chunk is Gaussian, the KL divergence will be zero. The KL divergence is calculated by:

Discussion
Constraining net rotation of the lithosphere is an important consideration for building absolute plate motion models. There have been several approaches to constraining LNR for the Earth. Most plate motion models produce some LNR. The LNR in these models arises because the modeled plate motions cannot fit geological and palaeomagnetic data with no rotation. Some studies choose to explicitly minimize LNR when fitting data in their plate motion models, for example, Tetley et al. (2019), although this assumption is not constrained by geodynamic considerations, primarily because geodynamic constraints on LNR were not previously available. A key consideration to determining the magnitude and direction of LNR in plate motion models is the choice of mantle reference frame used. Hotspots produced by mantle plumes are often used as a reference frame, with the assumption that the lithosphere moves with respect to quasi-stable plumes in the mantle. The choice of which plumes to include in the reference frame can alter estimates of net rotation (e.g., Gripp and Gordon (2002) and Doubrovine et al. (2012)). Because we have no reliable record of mantle structure in the geological past and the mantle convection simulations used to study convection evolution often use plate motion models as surface boundary conditions, it is difficult to assess how truly stable competing reference frames are.
Because we are studying convection models in this work, we can sidestep the problem of unstable reference frames. Lateral motion in StagYY is calculated relative to a no net surface rotation reference frame. It is then straightforward to integrate rotation in the mantle to give us LNR in a no mean mantle net rotation reference frame. This is impossible to calculate for the Earth, even at the present day because we have no way of reliably measuring lateral flow for the whole mantle. Perfect knowledge of the modeled LNR allows us to study the variations caused by both convection parameters and the evolution of mantle structures.
Previous studies have identified that lateral viscosity variations are required in geodynamic models to generate lithospheric net rotation, whether introduced through strain and temperature-dependent viscosity laws or by the addition of highly viscous continents (Hager & O'Connell, 1979;O'Connell et al., 1991;Bercovici et al., 2000). In our models, we find that the only first-order controls on LNR are choices made in the viscosity law, through the activation energy, activation volume, and presence of a viscosity jump at the 660 km transition zone. In agreement with previous studies, for example, Becker (2006), changes to viscosity in the upper mantle have a much greater effect on LNR, as seen in (Figure 5) where increasing a E E strengthens a positive correlation between upper mantle viscosity and LNR. As in Becker (2006) we posit that increased lateral viscosity variations promote surface mobility, allow more lateral motion in the asthenosphere and upper mantle, and potentially weaken the connection between lithosphere and mantle, promoting rotation. The correlation is stronger under the oceans, where more mobile lithosphere leads to a thinner thermal boundary layer and lower viscosity in the upper mantle, allowing more motion between the upper and lower mantle, and therefore more LNR in our models. This can be seen in the positive correlation between the near-surface temperature and LNR, where a hotter surface is the result of a more mobile surface. We find that the presence of continents has minimal effect on LNR, regardless of how thick they are, unless  3 a E E . This is broadly in agreement with previous studies, which either found continents to have a minimal effect (Gérault et al., 2012), or that continent are important where lateral viscosity variations are low (Becker, 2006;Zhong, 2001). In our models, the continents are passive and do not deform. They reduce the oceanic surface area available for forming ridges and subduction zones by 30% but otherwise do not significantly change the style of oceanic tectonics. In the previous studies where continents were found to make a difference, the viscosity contrast in the mantle was limited to around two orders of magnitude. The continents, therefore, provided the largest viscosity contrast, and slabs were not particularly viscous. In our simulations, at low values of a E E , the presence or absence of continents make more of a difference, in agreement with the comparable models in the previous studies. However, as soon as a E E is increased, the viscosity contrast of the stiff slabs in the mantle has a far greater effect on LNR than the presence or absence of continents. In line with this, we see a weak correlation between increased mid-mantle viscosity and magnitude of LNR, indicating that periods with more subducting slabs often have higher LNR. However, the correlation is not strong enough to claim that more slabs in the mantle always lead to greater LNR.
Our models are not stochastic, but they are highly non-linear. LNR is an emergent property of this non-linear system, where small changes in temperature or rheology at any point in their history can lead to global changes in the mantle structure via many non-linear interactions. The LNR, therefore, does not arise at random, but the interactions between different parts of the mantle flow which cause LNR are too complicated to be unraveled. We cannot identify any tectonic setting that reliably produces high or low LNR from studying the thermal and viscous structure of the mantle and lithosphere. We see that LNR operates on very different timescales to other parameters. LNR varies much faster than mean temperature and viscosity and does not change in line with changes to surface mobility.
With the relatively small number of models we have for this study, the complexity means that we cannot predict or invert for LNR at any particular point in time. With a much larger number of simulations (probably tens of thousands more), we may be able to use statistical methods, such as machine learning tools, to find characteristics of mantle flow that can determine the magnitude of LNR (Atkins et al., 2016;Gillooly et al., 2019). This is still not particularly practical for absolute plate motion models, because we have limited information on the tectonic arrangement of continents in the geological past. This is likely to limit the possibility of inverting for LNR.
However, despite the non-linearity at a statistical level, the distribution of net rotation is Gaussian, with the standard deviation dependent on the viscosity of the asthenosphere. The use of asthenospheric anisotropy to estimate the magnitude of LNR for the Earth (as in Conrad & Behn, 2010) therefore seems especially practical. LNR can be modeled as a Gaussian variable in absolute plate motion models or geodynamic simulations. This would allow LNR to vary if other data can be fitted better by using a higher value of LNR, with the knowledge that this is geodynamically reasonable.
We find that LNR is correlated over short time scales, up to about 5 Myr. For recent geological history, we have palaeomagnetic data with sampling rates greater than this, and we should therefore assume that absolute plate motion models using data with resolution below 5 Myr should have smoothly varying LNR. Over longer periods, this means that the jumps seen in models of LNR are not unreasonable and are therefore not necessarily a cause for concern (for a summary of the variability of LNR histories, see Figure 6 in Tetley et al., 2019). We also find that LNR is expected to have an approximately Gaussian distribution with sample lengths as short as 50 Myr.
We performed this study using a 2D spherical annulus geometry. The net rotation is therefore equivalent to that in the plane of a great circle. The velocities of convective motion in a spherical annulus case are generally similar to those produced by a full 3D simulation (Hernlund & Tackley, 2008), therefore we expect the LNR to have similar values. However, we cannot rule out that motion perpendicular to the plane may change the magnitude and mechanisms of net rotation, as the mantle flows around fully 3D slabs and continents. The range of root mean squared values of LNR in our models are generally somewhat lower than those proposed by absolute plate motion studies because we use lower than Earth-like viscosity variation, apart for models with high a E E ( a E E = 8-12), where we get values comparable to Earth (e.g., simulations 61-63, LNR = 0.18-0.289°/Myr, see Appendix A). That we succeed in producing Earth-like LNR in 2D model when we use high a E E suggests that we adequately model the key features responsible for producing LNR. Having large viscosity contrasts is also clearly a requirement for producing Earth-like values of LNR. Naturally, because our study was conducted in 2D, we cannot make any inferences about the likely stability of the pole of net rotation.
An interesting feature emerges when comparing the mean and the RMS values of LNR over the runtime of the models, presented in Appendix A. The models generally have mean LNR several orders of magnitude smaller than the RMS. This is expected because the LNR tends to oscillate around 0. This observation, combined with the observation that the correlation length of the models is not much more than 5 Myr may explain some of the mismatches between LNR predicted by absolute plate motion models and present-day estimates. The model of Torsvik et al. (2010) uses palaeomagnetic data for the last 5 Myr to produce their LNR estimate. Based on our models, this is long enough to begin to average over a large variation in LNR. In contrast, a single present-day estimate would be a sample from a Gaussian distribution with a standard deviation similar to the RMS values. It is, therefore, possible that both the larger, present-day values estimates in studies such as Gripp and Gordon (2002) and the smaller, time-averaged values from plate reconstructions can both be correct. Studies of seismic anisotropy are more likely to produce values closer to the RMS because shearing pre-deformed olivine in the mantle does not return it to an isotropic state. It, therefore, retains a record of the previous deformation even if the magnitude or direction of LNR changes.

Conclusion
For the first time, we present lithospheric net rotation measurements from fully self-consistent convection models with lateral viscosity variations comparable to those expected on Earth. We find that our models can achieve values in the same range as those expected for Earth, although only when the viscosity of the mantle is very temperature sensitive. We, therefore, conclude that the range of values previously published for Earth is geodynamically reasonable. The non-linearity of mantle convection makes predicting the magnitude and evolution of LNR very difficult, although we do find that increasing the number of slabs in the mid-mantle is likely, but not guaranteed, to increase LNR. Whilst we cannot predict the evolution of LNR, we do find that it is correlated over short time scales (up to 5 Myr). Beyond this, the distribution of values of LNR is Gaussian, especially if time scales longer than 50 Myr are considered. This provides a new geodynamic constraint on the range of amplitudes and correlation periods expected for Earth, which can be used to constrain future absolute plate motion models.

Data Availability Statement
A summary of the model data presented in this work is available to download in hdf5 format from Mendeley Data, https://doi.org/10.17632/z3g3zpktsk.1 (Atkins & Coltice, 2021), and includes model parameters, and the time evolution of net rotation, velocity, and viscosity.