Non-Linear Effects of Gravity Change on Mantle Dynamics

Abstract

Starting from observed variations of the length of day in the past 2.5 billion years, we calculate the corresponding gravity variation and evaluate the possible effects that such variation would have induced on the lithosphere and on mantle convection. The lithospheric stress induced by the observed gravity increase at the equator in 500 Ma is of the same order as the one associated with a convective cell. We model a gravity increase that would result in 1.3 × 10⁻² µGal or 1.3 × 10⁻¹⁰ m s⁻² over 10 years, a value that is not far from the detection thresholds of modern gravimeters. Gravity increase also affects mantle dynamics by improving convection efficiency. Our models demonstrate that gravity variations contribute to this phenomenon through faster and wider mixing. The results support a view where a convective system being subject to gravity changes responds through non-linear adjustments of its effective mixing at all scales. These findings contribute to a deeper understanding of how gravitational fluctuations may shape the behavior of Earth’s dynamic systems over geological timescales.

1. Introduction

Gravitational acceleration (g) is ubiquitously considered an invariant in all the geological and geophysical processes acting on the solid Earth, regardless of their spatial or temporal scales. Since it directly affects, among other parameters, the momentum conservation, Rayleigh number, and adiabatic heating of the mantle, gravity has direct effects on the very fundaments of mantle dynamics, regardless of their scale and dimensions [1,2,3,4], and affects all geological and geophysical processes from plate motion and tectonics [5] to crystal growth [6]. Gravity change in time is not included in the governing equations for geodynamical or geological processes since these processes act on different timescales than planetary gravity variation. Nevertheless, there is evidence suggesting that during the geological evolution of the Earth, gravity has changed. Mostly driven by tidal friction in the Earth–Moon system, tidal despinning, i.e., the decrease in Earth’s angular velocity, has affected the length of day (LOD) (Figure 1 in [7]), gravity [8], and more in general, the energy budget of the solid Earth [9,10], at least during the last 2.5 Ga (Figure 1a).

Despite the evidence for LOD variations (Figure 1 in [7]), retrieving the temporal trend of gravity at a given latitude (λ) through the past billions of years is not an easy task. This uncertainty is intrinsically connected to the geologic heterogeneity, dynamics and complex shape of planet Earth. While the geologic heterogeneity implies local variations of the resulting products and morphologies whose emplacement is strictly related also to gravity (e.g., the lateral extent and the compaction of a sedimentary deposit), the shape of planet Earth is responsible for global-scale variations of gravity, regardless of the geological process (e.g., latitudinal variations of g). We cannot take into account all possible combinations pertaining to the first uncertainty source, but we can simplify the shape of the Earth and estimate the latitudinal gravity variation rate (Δg(λ)) for the time interval covered by LOD data (see Figure 1).

Such simplification of the shape of the Earth can be numerically represented by Equation (1), which is similar to the latitudinal variation of gravity over a Clairaut spheroid [7] (Equation (10) in Appendix A), but allows us to account for temporal changes of m and f retrieved from LOD data (see Appendix A and Supplementary Figure S1).

$$g\left(\lambda , t\right)\approx {\displaystyle \frac{G M}{R^2}} \left[1+{\displaystyle \frac{1}{3}} f\left(t\right)-{\displaystyle \frac{3}{2}} m\left(t\right)+\left({\displaystyle \frac{5}{2}} m\left(t\right)-f\left(t\right)\right){sin}^2\lambda \right]$$

(1)

Equation (1) can be used to estimate “paleo-gravity” in the past at time t₁ using the adequate m and f retrieved from LOD data, and it can also be used to estimate the present-day g(λ, t₂) using the appropriate m and f (

Table A1. Parameters used for the estimates of g(λ, t₂) and Δg(λ). All values are taken from [13].

G (m3 kg−1 s−2)	6.67430 × 10−11
M (kg)	5.972300 × 10+24
R (m)	6.371000 × 10+6
present-day m	3.46775 × 10−3
present-day f	3.35281 × 10−3

in Appendix A). Finally, we can compute the gravity variation rate at each latitude Δg(λ) shown in Figure 1b:

$$\mathsf{\Delta }g\left(\lambda \right)=g\left(\lambda , t_2\right)-g\left(\lambda , t_1\right)$$

(2)

Δg(λ) represents the average gravity variation rate in the past 2.5 Ga accounting for angular velocity and flattening changes due to Earth’s tidal despinning. Figure 1b shows that Δg > 0 in the +62/−62 latitude range, with maximum average values equal to +1.30 mGal Ma⁻¹ (+1.3 × 10⁻⁵ m s⁻² Ma⁻¹) found at the equator, while at higher and lower latitudes we found smaller or negative Δg. Such gravity variation approximately corresponds to ~4 m height loss in one million years. As a comparison, geoid adjustment at the mid-Atlantic ridge is ~0.16 m Ma⁻¹ [11].

These values should be considered as the gravity variation due to changes of Earth’s flattening and angular velocity in the last 2.5 Ga, assuming that these are the only contributors to the observed LOD changes—i.e., excluding Earth’s mass and average radius variations. Furthermore, we exclude contributions from tidal correction factors which, at the timescale we are considering (Ga), may be regarded as transients [12]. We also assume other external or internal contributors to variations of Earth’s moment of inertia as transients since these can be regarded as temporally (e.g., glaciation) or spatially (e.g., subduction) transient if compared to the temporal (hundreds of Ma) or spatial (planetary scale) scales we are modeling.

The combination of flattening variation and angular velocity decrease [13], by reflecting a tidally driven slowing and change in the shape of the planet, produces opposed effects between equatorial and polar latitudes, in agreement with relative weights of flattening and angular velocity on Δg(λ) (see Equation (1)). Effects related to the latitudinal variations of gravity are represented here as instantaneously generated by angular velocity decrease, assuming a planet without internal friction. This is obviously not the case, but considering that the gravity increase at the equator and decrease at the poles are one-directional, since no alternating effects are expected (like it would be in the case of Earth body tides), any temporal delay between the decrease of angular velocity, the corresponding LOD increase, and the corresponding adjustment of Earth’s ellipticity and the following gravity variation would not affect the magnitude or the position of the Δg estimated here since, whatever the magnitude of such delay, it would equally affect the entire time series we consider.

The LOD data we use to estimate m and f temporal trends show steeper diminishing rates and denser data points between 0 and ~0.5 Ga before the present than in the 0.5–2.5 Ga timespan (Figure 1 [7,9,14,15]). In turn, this implies that f and m decreased more rapidly in the last 500 Ma than in the previous 2 Ga, and that the estimates of f and m in the last 500 Ma are more reliable than those in the previous 2 Ga. Hence, the equatorial gravity increase rate of 1.3 × 10⁻⁵ m s⁻² Ma⁻¹ that we calculated over the entire LOD data series (2.5 Ga) represents a lower bound if the 0–0.5 Ga period is considered, whilst it represents an upper bound in the 0.5–2.5 Ga period. As an example, the Permian-to-recent equatorial gravity increase was ~2 mGal Ma⁻¹ [7]. However, since focusing our models on a particular epoch would require some other constraints to be included in the model (e.g., geometrical and/or petrological), considering that we want to keep our models as simple as possible and free from any temporal constraints, in the following we use the average 1.3 × 10⁻⁵ m s⁻² Ma⁻¹ equatorial gravity increase that we calculated over 2.5 Ga.

2. Geodynamic Implications

The questions arising from the Δg estimates in Figure 1b are the following: What are (if any) the geodynamical consequences of such changes in gravity? And are such effects observable and/or quantifiable?

Working on the actually observable geodynamic domain, some authors proposed some relations between LOD variations and seismicity distribution [10,16] or with GPS measurements [17]. But what was the Δg effect on the Earth 2.5 Ga ago? Given the lack of evidence, addressing this question requires several assumptions. In the following, we evaluate the effects that a positive gravity rate, such as the one found over equatorial latitudes, would have on lithospheric and mantle dynamics.

A 200.5 km thick continental lithospheric column (cratonic) at the equator, 2.5 Ga ago, is taken as consisting of 160 km of mantle (with density ρ_m = 3300 kg m⁻³) and 40.5 km of crust (with density ρ_c = 2800 kg m⁻³) including a minimal topography height (h₀) of 0.5 km. If we consider a present-day equatorial gravity of 9.78 m s⁻² and the average gravity increase at the equator from Figure 1b (1.3 × 10⁻⁵ m s⁻² Ma⁻¹), we can estimate a paleo-gravity of 9.75 m s⁻². Such lithospheric structure would result in a gravitational potential energy (GPE) per unit area (Γ) of 611 TN m⁻¹, representing the potential energy of the lithospheric column. This, in turn, would produce a vertical depth-averaged stress (σ_zz) acting on the underlying mantle of 3.048 GPa (see Appendix A and Supplementary Figure S2).

The balance between GPE and normal deforming stresses results in the observed deformation of the continental lithosphere [18,19,20,21]. Local variations of GPE have often been investigated to link surface evidence (e.g., geoid undulations, GPS velocity maps) to stresses transmitted through the lithosphere from deeper sources (e.g., convective or tectonic sources) [20,21,22].

However, in the case we are considering—i.e., a stagnant lithospheric lid not disturbed by convection or tectonic perturbations—the only source for variations of GPE, and in turn, in σ_zz, is g. If we fast-forward in time the equatorial lithospheric column by 500 Ma, without external perturbations, the only parameter affected by the elapsed time is g, which now is 9.7565 m s⁻² with a Δg of 6.5 × 10⁻³ m s⁻². If we assume the density values and the geometry of the lithospheric column as constants, the observed Δg results in a vertical stress increase Δσ_zz of ~2 MPa (i.e., a Δσ_zz rate of ~4 kPa Ma⁻¹). Eventual contributions related to density variation can be assumed negligible (see Appendix A).

The resulting vertical stress increase is uniquely related to gravity increase and is directly acting on the underlying mantle, and its magnitude is comparable with the stress associated with a convective cell of radius 10⁶ m, velocity of 10⁻⁹ m s⁻¹ and average viscosity of 10²¹ Pa s [23].

3. Computational Mantle Convection Models

In the following, through convective numerical models and entropy estimates of the mixing process, we estimate effects of gravity variations on the mantle. The aim of this work is to evaluate if and to what extent a gravity variation could affect mantle dynamics, which are extremely slow. Accordingly, in the following, rather than moving in space (from equator to poles) in a single temporal snapshot, we will attempt to model dg/dt contributions over hundreds of Ma at the equator, where the gravity increase in time was more pronounced (Figure 1b).

If we move from the unidimensional lithospheric case towards a wider two-dimensional model of the lithosphere and the mantle, the effects of Δg can be evaluated in a more compelling view. Figure 2 shows the 2D visco-elasto-plastic model with 100 km of lithosphere overlying 2900 km of mantle. The modeling assumes a laterally isothermal starting state (i.e., no initial lateral temperature variations, see Supplementary Figure S3).

The visco-elasto-plastic models shown in Figure 2 were computed using the Mantle_convection.m code presented and discussed in [3]. In agreement with the original code, the models shown here were created within a 3000 km × 3000 km domain (51 × 51 nodal points and 40,000 randomly distributed markers). Free-slip boundary conditions are prescribed on all boundaries. No heat flow is allowed along the vertical boundaries, while constant temperatures are set at the top (273 K) and at the bottom (3500 K) boundaries, 0 and 3000 km depth, respectively. The input initial thermal state was set in order to avoid lateral temperature changes at modeling start (see Supplementary Figure S3). Such initial setup was chosen in order to avoid any possible contribution to the modeled convection from pre-set lateral thermal variations. This may sound unrealistic for a real case convective mantle, but we preferred to remove any possible contributions from other sources to focus our investigation on dg/dt effects on the convective system. We discuss effects of lateral thermal variations in Section 6. Furthermore, the initial gravity gy was set to 9.75 to simulate a model starting 2.5 Ga ago. Finally, to model gravity increase, we also include a gravity increase function to reflect the observed equatorial gravity increase in time. Assuming that we are modeling an E-W section over the equator, the original code was modified to include the gravity increase rate (dg/dt) of 1.3 × 10⁻⁵ m s⁻² Ma⁻¹ retrieved from LOD variations (Figure 1b).

The modeling is based on the second invariants of the deviatoric strain rate tensor ($\dot{{\epsilon }_{II}}$) and deviatoric stress (σ_II) of the dry olivine non-linear flow law (Table 21.1 in [3]),

$$\dot{{\epsilon }_{II}}=A_{D }{{\sigma }_{II}}^n exp \left(-{\displaystyle \frac{E_a}{R T}}\right)$$

(3)

where, according to [3], the material parameter A_D is 2.5 × 10⁴ M Pa^-n s⁻¹, the stress exponent n is 3.5 and the activation energy E_a is 532 kJ mol⁻¹; R is the gas constant and T is temperature.

Thus, the effective viscosity η follows the power law [23]

$$\eta ={\displaystyle \frac{{\sigma }_{II}}{\left(2 \dot{{\epsilon }_{II}}\right)}}={\displaystyle \frac{{\sigma }_{II}^{(1-n)}}{\left(2 A_D\right)}} exp \left({\displaystyle \frac{E_a}{R T}}\right)$$

(4)

In all models, mantle viscosities range between 10²⁰ and 10²³ Pa s, either if a gravity rate or constant gravity value is used.

Following the original setup of the modeling code, we do not consider gravity variations with depth in the mantle. Although this may look like an over-approximation, it does not affect the results of our modeling since the predicted gravity increase within the mantle (e.g., the PREM model [24]) would equally affect both the model with and without the gravity rate we model (Figure 2), increasing the baseline gravity in the lower mantle, without any effect on the following entropy estimates (Figure 3b,c).

Despite the similarities of the modes regarding the general trend of mantle mixing, some slight differences are observable between the model, including a gravity change in time and the coeval model without gravity increase. From the mixing observed in both models in the upper mantle, down-going thermal plumes generated earlier in the model including gravity increase (white dashed squares in Figure 2). On the contrary, up-going plumes are not affected by the relatively higher contrasting gravity. Such differences are clearly observable up to ~170 Ma after modeling started, when mixing becomes more intense in both models. In the long-period modeling, up to ~600 Ma, when gravity increase is included, the propagation of mixing towards intermediate depths and to the core–mantle boundary is more effective both qualitatively (white dotted squares in Figure 2) and quantitatively, as we will show in Section 4, where the geological entropy (also known as configurational entropy according to [25]) of the convective models will be estimated from the mixing of blue and violet particles (Figure 3).

In general, even in simple models such as those shown in Figure 2, including the gravity variation results in faster and more effective downwards propagation of thermal anomalies. From the modeling perspective, this is unsurprising since the gravity directly affects both the Stokes and the temperature equations, solved through a free energy minimization approach using a visco-elasto-plastic modeling code [3]. The reader is referred to [3] for a complete description of the code. In the next section, we will use these models to assess the eventual effects of the gravity change per million years (1.3 × 10⁻⁵ m s⁻²) on the mixing efficiency of the convective system.

4. Estimating the Geological Entropy of the Mantle

To quantitatively evaluate the significance of the changes induced by the gravity variation, we estimate the geological entropy of the convective models using directional entrograms [26,27]. Geological entropy measures the disorder of a system through the quantification of the uncertainty that a random source of information contains—i.e., the Shannon information entropy (H) [28]. For completely ordered systems, H → 0, while for chaotic systems, H → 1. It should be noted that the Shannon entropy is mathematically analogous to the Boltzmann–Gibbs entropy and thus can be used to describe the thermodynamic state of a system.

The Shannon information entropy H of a system M containing N categories (or events) m_i (i = 1, …, N) is given by

$$H\left(M\right)=-\sum _{i=1}^{N}P\left(m_i\right)\left[-\ln P\left(m_i\right)\right]$$

(5)

where P(m_i) is the probability of the ith category (or event) to occur [28].

Directional entrograms shown in Figure 3 are calculated from the local entropy H_L as defined in [27,29]. In their algorithm, the entropy is calculated in randomly placed local subdomains of the system (searching boxes) that in the 3D system have dimension l = l_x i + l_y j + l_z k. The size of l is defined by the unit vectors i = (1, 0, 0), j = (0, 1, 0) and k = (0, 0, 1) and the scalar components l_x, l_y and l_z, while l_z = 1 in 2D systems. The smallest searching box in 2D systems is a square (Figure 3) where l_x = l_y = 1 times the grid size of the system (10 km in our convection models).

The local entropy H_L(l, j) of the searching box of size l (Equation (6)) is calculated N_MC times (j = 1, …, N_MC).

$$H_L\left(l, j\right)=-\sum _{i=1}^N{p}_{L, i}\left(l, j\right)\left[\ln p_{L, i}\left(l, j\right)\right]$$

(6)

where p_L,I represents the proportions of the categories within the searching box (e.g., the fraction of cells with blue or violet color in the left panels of Figure 3). The N_MC value represents the number of times that the searching box changes location, maintaining its size, every time starting from a random point of the system, and calculates the local entropy. Once N_MC calculations are achieved, the iteration restarts with a larger searching box. Following [27] in our calculations, we always maintain N_MC ≥ 500.

At each calculation step, H_L(l, j) is normalized over the entropy H₀ of the entire system (0.69 in our case) to obtain a local relative entropy H′_R (l, j) = H_L (l, j)/H₀. Finally, an average relative entropy of the searching box H_R(l) is calculated:

$$H_R\left(l\right)={\displaystyle \frac{1}{N_{MC}}} \sum _{j=1}^{N_{MC}}{H}_R^{\prime }(l, j)$$

(7)

In a chaotic system, H_R(l) → 1 since H_L (l, j) → H₀.

The entrograms shown in Figure 3 plot the N_MC-averaged local entropy H_L(l) and the relative entropy H_R(l) against the size of the searching box (l) expressed in the number of cells. In a 2D system such as the one modeled in this work, changing the i and j unit vectors of l results in tuning the horizontality and verticality of the searching box, and, in turn, allows us to calculate the local and relative entropies along the horizontal (red plots in Figure 3) and vertical (green plots in Figure 3) directions. To allow the searching box to remain within the analyzed volume, its maximum extent is half of the maximum number of cells in the searching direction—i.e., 150 cells or 1500 km both in the horizontal and vertical direction. The minimum extent of the searching box, i.e., when l_x = l_y = 1, gives the entropy of the smallest element of the system (H_L0).

Since the relative entropy H_R(l) of the searching box is representative of the spatial disorder along the measured horizontal or vertical direction, the normalized entropic scale nH_SX [27]

$$n{H}_{SX}={\displaystyle \frac{1}{\max l_X}} {\int }_0^{\infty }\left[1-H_R\left(l_X\right)\right]dl$$

(8)

represents the spatial order along the horizontal direction.

If we consider the convection models shown in Figure 2 as “snapshots” of an evolving convective system, by measuring the geological entropy at each step, we can measure the efficiency of the convective system—i.e., its ability to distribute chaos over the initial order or, in other words, its ability to mix the initially layered particles.

In a system such as the one modeled in Figure 2, where the thermodynamic equilibrium is maintained through Gibbs free energy minimization [3], we expect the general entropy (H₀) to remain constant throughout the modeling phases. However, the entropy of small samples of this volume and the entropy along sampling directions should increase according to convective mixing. In the initial setup, when the horizontal layering of the modeled volume is undisturbed (upper central image in Figure 2), we measure directional entropy along the horizontal direction (red plot in Figure 3a), which is significantly lower than the directional entropy along the vertical direction (green plot in Figure 3a). At a later modeling step, the entropy along the horizontal direction will increase due to vertical and lateral flow—i.e., because of the convective displacement of the modeled particles. By tuning the size of the sampling box, the entropy can also be calculated for small portions of the system, providing a local entropy estimate (H_L). Finally, the entropy of the smallest element included in the convective system (H_L0) can be estimated by reducing the sampling box for H_L calculation to the size of the cell of the modeled system (10 km × 10 km).

When the directional entropies are equal, i.e., when the directional entrograms overlap, the system is isotropic along those directions. In general, the increasing disorder of the system is reflected in the increasing entropy of its smallest cell, towards values comparable to the global entropy of the system (H_L0 → H₀).

5. Results

The entrograms shown in Figure 3 depict the evolution of the local entropy H_L in the vertical and horizontal directions as mixing increases over time. Entropy along the horizontal direction before convective mixing is significantly lower (coherence is higher) than the one measured along the vertical direction, regardless of the size of the sampling box. This is due to the horizontal layering imposed in the initial setup. Note that such layering is not representative of any compositional or thermal variations but is included in the modeling to visualize the convective distribution in the following temporal steps. As mixing goes on in time (Figure 3b), the entrogram along the horizontal direction (red curve in the plot) tends to match the vertical one (green curve in the plot). As time goes on, we can see that the distance (size of the sampling box) required for the local entropies to reach the system entropy (H₀) becomes smaller; i.e., the entropy increase spreads to smaller volumes of the system. After a significant amount of time (Figure 3c), the directional entrograms become more and more similar because the mixing has almost reached the maximum disorder at all scales (H_L0 → H₀).

All of the above is not surprising since we are dealing with a system that is being convectively mixed. However, if we include a gravity increase rate in time (continuous lines in plots of Figure 3b), and compare these entrograms with those calculated in a convective system without gravity changes (dashed lines in plots of Figure 3b), we see that the former provides slightly higher entropy values both on the horizontal and vertical entrogram (Figure 3b). While these differences are somehow observable in the 200 Ma time step, after 600 Ma of convective mixing, the directional entrograms become almost indistinguishable (Figure 3c).

We interpret these differences as representative of contributions to the efficiency of the convective mixing provided by the gravity increase. Such contributions allow for a slightly faster propagation of entropy increase to the smaller components of the convective system. This interpretation is supported by the entropy values of the smallest component of the modeled system—i.e., the 10 km × 10 km cell (H_L0). This parameter can be regarded as an index of propagation of the entropy increase due to mixing. The faster the entropy of the smallest element of the model increases, the more efficient the mixing is. We can see that accounting for gravity increase in the convective model after 200 Ma provides H_L0 = 0.37, while neglecting the gravity increase provides an H_L0 = 0.35.

Since the geological entropy ranges between 0 and 1, we can estimate that accounting for gravity increase in the convective model provides some 2% increase of the entropy of the cell. In our view, such contribution is significant when compared with the 26% entropy increase observed in the 200 Ma model with constant gravity (H_L0 increases from 0.09 to 0.35, Figure 3b).

Such contribution is also observable if we consider the normalized entropic scales along the two sampling directions (nH_SX and nH_SY). Since these parameters measure the spatial order along the horizontal and vertical directions, the convective mixing lowers the initial order along the horizontal direction (nH_SX) towards zero. On the other hand, the initial lower order along the vertical direction (nH_SY = 0.04) is almost preserved throughout the modeling if constant gravity is maintained (nH_SY = 0.03), but still, some contribution can be observed if dg/dt is accounted for (nH_SY = 0.02). Also, in this case, accounting for gravity changes in the convective model results in lower values of nH_SX in respect to models where the gravity change is not included.

All computed entrograms are grouped in one single plot in Figure 4. These plots are computed to map the spatial and temporal evolution of H_L along the horizontal direction (color-coded values) against the time since the convective model started (horizontal axis) and the size of the searching box (vertical axis). Left plots account for entrograms computed considering dg/dt, while right plots show entrograms computed with constant gravity.

The color variation parallel to the abscissa depicts the temporal evolution of the Shannon entropy of a given searching box along the horizontal direction, while a color variation parallel to the ordinate depicts the spatial distribution of the Shannon entropy at a given epoch after the convective model has started. For example, the plot shows that in a searching box of 1000 km size, H_L reaches 0.65~400 Ma after convection starts if gravity is kept constant, while it reaches the same Shannon entropy ~100 Ma earlier if gravity increase is considered (Figure 4a,b). Furthermore, the minimum size of a searching box with the highest entropy is somehow representative of the lateral spreading of the entropy increase due to mixing: the smaller the size of this searching box, the more widespread the entropy increase. For example, at 500 Ma, the Shannon entropy is ≥0.65 for all searching boxes > 440 km if gravity increase is not considered (Figure 4d), but this spatial threshold decreases by ~100 km (entropy increase is more widespread) if the gravity increase is accounted for (Figure 4c).

Regardless of the plots we analyze in Figure 4a or Figure 4b—i.e., either if we consider a gravity increase trend in time or not—we note that within 100 Ma after convection begins, the entropy increase is widespread, with H_L ≥ 0.5 for all searching boxes larger than 100 km. This observation depicts a fast entropy increase in the early stages of convective mixing where the contribution of gravity increase is undetectable. However, the more the convective mixing goes on, the more the contribution due to gravity increase is revealed since the propagation of isentropic contours in the plots with constant gravity values is slower. Considering the isentropic contours of 0.60 and 0.65, i.e., the H_L values that are closer to 90% H₀ (0.62), their spatial and temporal distributions are noticeably different: when the gravity increase is accounted for, it results in earlier entropy increase at all the observed scales (Figure 4).

Bearing in mind the general significance of entropy, and considering that the Shannon entropy has been used in a wide range of applications to evaluate the level of uncertainty of a system’s random variable [30,31,32,33,34], in the theoretical setup shown here, we can consider the changeof entropy of the smallest element of the convective system (H_L0) as representative of the efficiency of the system. The more the convection is active, the more mixing it will achieve, and the higher the entropy of a small sample of the convective volume will be. Such view holds both when considering a single time interval (the higher the entropy of the cell, the higher the mixing efficiency of the system) and if a temporal trend of the entropy variation of the cell is considered (the faster its entropy rises, the more effective the mixing is). In this view, the time rate of H_L0 can represent the evolution of the convective mixing from an initial state where the order is higher (H_L0 → 0, nH_SX → 1) towards a later stage where the convective flow has increased the entropy at all scales towards the general entropy of the system (H_L0 → H₀, nH_SX → 0).

Figure 5 shows the temporal trend of the entropy (H_L0) and of the order (nH_SX) along the horizontal direction of the smallest element of the convective model—i.e., the 10 km × 10 km cell, between 0 and 608 Ma after convection modeling begins. Figure 5a shows that ~120 Ma after the model started, the data series including the gravity increase (black dots) diverges from the data series not accounting for it (red squares). Such divergence is maximum at ~300 Ma, where the H_L0 value accounting for dg/dt is 0.428 while H_L0 with constant gravity is 0.407. The divergence is almost null at ~608 Ma, where the difference between H_L0 with and without dg/dt is 0.001. A similar trend is observed for the nH_SX plot in Figure 5b, with maximum gaps between the two models of 0.01 found between ~100 and ~400 Ma.

In summary, the gravity increase estimated from LOD data provides adjunct vertical loads at the base of the lithosphere comparable to those associated with a convective cell. Furthermore, accounting for gravity variation in time provides some contribution to mixing efficiency of a convecting volume. In particular, when accounting for gravity changes, we found a faster propagation of the Shannon entropy increase within the convective model and a corresponding faster decrease of order (Figure 4 and Figure 5).

We are aware of the limited representativeness of the simple models proposed here, and the purpose of this work is not to represent a real planetary-scale convective mantle. However, in the following discussion, we aim to generate some interest in the possible effects of gravity variations, since the contributions related to gravity increase are quantifiable and may be significant in the long-term (>100 Ma) evolution of the convective system of the Earth or other rocky planets.

6. Discussion and Conclusions

Bearing in mind that our models start from an equally abundant distribution of blue and violet particles (Figure 3a), recalling that the Shannon entropy is mathematically analogous to the Boltzmann–Gibbs thermodynamic entropy, and that the gravity increase rate is kept constant throughout the modeled time, it follows that the time required for the Shannon entropy of the smallest cell to increase is governed by the convective mixing: the more efficient the mixing, the less time it will take to increase the entropy throughout the system. In this way, we can regard the convection as a disordering phenomenon that from a theoretical initial time before convection (i.e., when blue and violet particles are ordered, upper central image in Figure 2) distributes entropy to the smallest portions of the system towards an isentropic state—i.e., H_L0 → H₀. If all the above holds, including gravity increase in the evolution of the convective system allows for faster spreading of the entropy perturbation at all the observed scales, from the smallest 10 km × 10 km cell to the largest (Figure 4 and Figure 5).

It has been estimated that to reach effective mixing, a single-layer convective mantle requires ~500 Ma, with the convective mixing being responsible for the adiabatic (isentropic) condition in ~90% of the mantle (page 349 in [35]). Since no thermodynamic system is perfectly insulated, we can assume that equilibrium has been reached when H_L ≥ 0.9 H₀—i.e., when 90% of the general maximum entropy of the system (H₀) has been spread through convective mixing. Given the link between convection, entropy, mixing and adiabaticity of the mantle, and considering that in a chaotic dynamic system such as the convecting mantle [36,37,38] even a small contribution to entropy increase may be significant at all scales, we conclude that gravity increase allows for faster and wider effective mixing and directly affects the thermodynamic evolution of the convective system towards adiabatic conditions (Figure 4). Consequently, we suggest that a convective system being subject to gravity variations responds to such changes at all scales through non-linear adjustments of its effective mixing. In the following, we provide some arguments and tests to support this view.

The traditional way to describe the vigor of a convective system is the dimensionless Rayleigh number (Ra), which describes the ratio between buoyancy and diffusive forces acting on the convecting material (Equation (9)). The larger the Ra value, the more effective the convective heat transfer compared to diffusion. The Ra of a convective system of thickness D = 2.85 × 10⁶ m and average density ρ_m = 4000 kg m⁻³ heated from below is given by [11]

$$Ra={\displaystyle \frac{{\rho }_m^2 g \alpha \left(T_{bot}-T_{top}\right)D^{3}Cp}{\eta k}}$$

(9)

where g is gravity, α is the average thermal expansion coefficient equal to 1.6 × 10⁻⁵ K⁻¹ [23], T_bot (3700 K) and T_top 1500 K denote the temperatures at the base and on top of the convective system, respectively, Cp is the specific heat of the mantle equal to 1.3 × 10³ J kg⁻¹ K⁻¹ [23], η is the viscosity and k is the average thermal conductivity of the mantle equal to 8 w m⁻¹ K⁻¹ [23].

By maintaining constant the parameters quantified here (D, ρ_m, T_bot, T_top, α, Cp and k), we extract η from the convective models after 5 Ma and after 608 Ma of convection both with and without gravity increase. Viscosity values are extracted from the center of the modeled domain (i.e., from the 25th column of the 51 × 51 nodal points, see Section 3) and averaged, and the resulting value is used to calculate Ra. Vertical gravity g is 9.75 m s⁻² for calculations referring to constant-gravity models and 9.7579 m s⁻² after 608 Ma accounting for gravity increase.

After 5 Ma from model start, an average mantle viscosity at the center of our model equal to 6 × 10²² Pa s provides an Ra of 3.5 × 10⁵ no matter whether the gravity increase is considered or not. After 608 Ma of mantle convection with constant gravity, the average viscosity at the center of the model is 4.2 × 10²⁰ Pa s and the corresponding Ra is 4.9 × 10⁷. If the gravity increase is taken into account, after 608 Ma, the average viscosity at the center of the model is 3.2 × 10²⁰ Pa s and the corresponding Ra is 6.5 × 10⁸ with a ~30% increase of Ra due to gravity increase directly affecting both Ra (Equation (9)) and the viscosity of the resulting model. Such effect of the gravity increase on Ra provides further proof of the non-linear contribution of dg/dt to the mixing efficiency of the convective system.

In Figure S3, we compare the final (after 608 Ma) thermal state of the convective model and its viscosity values when dg/dt is taken into account (Figure S3b and Figure S3e, respectively), against the thermal state (Figure S3c) and the viscosity values (Figure S3f) when modeling is performed with constant gravity. The comparison provides similar viscosity and temperature values across both models, but when the velocity of the mantle flow is plotted (white arrows in Figure S3e,f), we notice relatively higher values (up to 40 mm y⁻¹) in the model accounting for dg/dt (Figure S3e). Such velocity values are in general agreement with those retrieved from slab sinking rates [39]. All the above suggests that the increase of mantle mixing efficiency (Figure 4 and Figure 5) is not caused by the modeled viscosity values.

Furthermore, in order to evaluate possible effects of the initial thermal state of the mantle on the entropy, we run the same convective models with an initial thermal model including some lateral perturbations (Supplementary Figure S4) and test this thermal model under both dg/dt or constant-gravity conditions. Although the effects of a different initial thermal state result in different viscosity values within the modeled domain (compare Supplementary Figures S3e and S4e), when the entropy plot is computed (Supplementary Figure S5), the contribution due to gravity is clear and coherent. This test demonstrates that whatever the uncertainties introduced by other factors (e.g., temperature) and their consequences, the gravity contribution on entropy spreading remains because, similarly to the model discussed above (Figure 4), the measured entropy highlights the positive contribution of dg/dt to its spread through space and time (Supplementary Figure S5).

Finally, to test the possibility that the entropy spreading within the convective system was due to some other uncontrolled “butterfly effect” [40,41,42] rather than being triggered by the small dg/dt perturbation introduced in the convective mixing, we run the very same model—i.e., the one including gravity increase—a second time with the same identical setup as the one discussed above (Figure 2, Figure 3, Figure 4 and Figure 5). In Supplementary Figure S6, we compare these two runs. Despite the entropy spreading within the first 150 Ma being slightly different between the models, the spatial and temporal trend of the maximum entropy (H_L = 0.65) converges to similar conditions between the two models, supporting a consistent contribution of dg/dt on entropy spreading over long time and distances. These tests (Figures S4–S6) demonstrate that the spatial and temporal evolution of the entropy of the two models, i.e., the one with constant gravity and the one with dg/dt, consistently diverges because of the small perturbation introduced by the gravity change.

The gravity variations discussed here are clearly not acting alone. Mantle dynamics can be affected by several parameters ranging from compositional to physical variations that are often combined with each other. To evaluate the relevance of each parameter in comparison with others is the main focus of 2D, 3D and 4D geodynamic modeling [1,3].

The increase of the precision and accuracy of gravimeters [43] may continue towards levels that allow gravity variation to be directly observed. Such a target (the gravity increase we model would result in 1.3 × 10⁻² µGal or 1.3 × 10⁻¹⁰ m s⁻² over 10 years) is not too far if we consider that long-term stability of absolute quantum gravimeters is <1 µGal [44]. In the future, some possible further contributions to dg/dt detection may come from planetary-scale 3D modeling and transformations [45,46] of high-accuracy (≈1 µGal) and high-precision satellite gravity observations. In such a case, the dg/dt contribution we model could represent a residual once lithospheric sources are removed.

The next investigative step regarding the effects of gravity increase on geodynamics should involve the understanding of any possible relation between the dg/dt rate and the forces driving plate motions. To understand such eventual relation(s) requires addressing the variable contributions due to gravity, beyond the common definitions of ridge push (RP) and slab pull (SP). There is a direct proportionality between g and RP [11] and between g and SP [5], so one would expect that increasing gravity would correspond to an increase of RP and SP, and in turn, to more efficient plate dynamics, similarly to that observed for the convective mantle (Figure 4). However, assessing whether the contribution related to the gravity increase in the +62/−62° latitudes has been significant in regard to plate motions would require a denser and possibly wider (older than 2.5 Ga) time series of data constraining the dg/dt rate estimates, mostly to understand its effects on compensation models (i.e., geoid adjustments and isostasy). Furthermore, some stronger temporal constraints would also be required regarding the beginning of plate tectonics. In a frame where the relative importance of the different forces driving plate motions is still a matter of debate [5], such an investigation would be beyond the scope of the present work.

Our current knowledge of mantle rheology and plate motion history carries uncertainties that would prevail over the dg/dt effects modeled here. The effects of gravity change are relatively minor and may be difficult to distinguish from those of other parameters. Nonetheless, the main point of this work is that they do exist and have physical consequences on mantle dynamics, even if they are, at present, difficult to detect.