Artificial Intelligence in Geosciences

Modern geodynamics is based on the study of a large set of models, with the variation of many parameters, whose analysis in the future will require Machine Learning to be analyzed. We introduce here for the first time how a formulation of the Lattice Boltzmann Method capable of modeling plate tectonics, with the introduction of plastic non-linear rheology, is able to reproduce the breaking of the upper boundary layer of the convecting mantle in plates. Numerical simulation of the earth’s mantle and lithospheric plates is a challenging task for traditional methods of numerical solution to partial differential equations (PDE’s) due to the need to model sharp and large viscosity contrasts, temperature dependent viscosity and highly nonlinear rheologies. Nonlinear rheologies such as plastic or dislocation creep are important in giving mantle convection a past history. We present a thermal Lattice Boltzmann Method (LBM) as an alternative to PDE-based solutions for simulating time-dependent mantle dynamics, and demonstrate that the LBM is capable of modeling an extremely nonlinear plastic rheology. This nonlinear rheology leads to the emergence plate tectonic like behavior and history from a two layer viscosity model. These results demonstrate that the LBM offers a means to study the effect of highly nonlinear rheologies on earth and exoplanet dynamics and evolution.


Introduction
The discovery of plate tectonics emerged from a huge amount of new data on seafloor bathymetry and a courageous new approach to analyze them.In the next decade, our way to use modeling to understand geodynamics will have to evolve by combining artificial intelligence with high-performance computing (HPC) simulations (Morra et al., 2021).While the number of geoscientists who use HPC is increasing, there is great need for fast and reliable geodynamic tools for producing the large number of training and verification samples necessary to train machine learning algorithms.New numerical methods for geodynamics must be carefully crafted to maintain their performance on parallel computers and be tested for non-linear rheologies which are responsible for the fragmentation of the surface in tectonic plates.
Geological processes such as mountain building, geomorphological evolution of river deltas, require a knowledge of the past, and hence its past history.Therefore, an important criterion for understanding the history of plate tectonics and the underlying mantle convection is a rheology which can preserve its history.Both viscoelastic and plastic rheology (Karato, 2008) can produce this attribute, with viscoelastic rheology limited to a time interval of at most tens of million * Corresponding author.
The Lattice Boltzmann Method (LBM) is an alternative approach to classical numerical solutions to PDE's for modeling fluid dynamics based on simulating the Boltzmann Equation on a discrete lattice.Since the first inception as the Lattice Gas Method in which it was proven that the Navier-Stokes equations are recovered at the macroscopic scale by modeling discrete particles moving and colliding on a discrete lattice (Frisch et al., 1986), the more efficient Lattice Boltzmann Method (LBM) was developed (Chen and Doolen, 1998) in which particle number densities are moved on a discrete lattice with collisions being achieved by relaxing the number densities towards the Boltzmann equilibrium distribution (Bhatnagar et al., 1954).The LBM, has since been used to model a wide range of fluid dynamics problems (Succi) including mantle convection (Mora andYuen, 2017, 2018a,b), magneto-hydrodynamics (Dellar, 2002), multiphase flow in porous media (Huang et al., 2015;Mora et al., 2021) and plastic flows (Falk et al., 2015).In the following, we overview the thermal Lattice Boltzmann Method which can be used to simulate mantle and exoplanet convection, and present 2D simulations in a two layer model with a highly nonlinear rheology and demonstrate that this rheology leads to the emergence of plate tectonic like behavior in the model.This is consistent with simulation results using classical methods (Tackley, 2000;Gerya, 2019) which require a nonlinear rheology such as plastic or dislocation creep to model convection with plate tectonics.However, we use a much more nonlinear rheology than the standard plastic or dislocation creep rheologies with an effective viscosity -above a yield stress -that is inversely proportional to strain rate to the fourth power (cf.plastic and dislocation creep rheologies respectively have effective viscosities inversely proportional to the strain rate to the power of one and to a power that is less than one).Use of this extremely nonlinear rheology demonstrates the LBM's potential to be used to study the effect of different highly nonlinear rheologies on earth and exoplanet mantle convection dynamics and evolution.Machine Learning (ML) and Deep Neural Networks have proven to be extremely useful scientific tools (Baldi, 2021).Previous work on machine learning applied to the geosciences has shown that models like the ones presented here can be organized as data to be analyzed using Convolutional Neural Networks (Dye and Morra, 2020) and other machine learning algorithms (Qingkai et al., 2021).

The thermal lattice Boltzmann method
The thermal LBM models two distributions   and   representing the mass density and energy density moving one lattice spacing along the  directions (i.e.orthogonally and diagonally) of a cartesian lattice.After these number densities move by one lattice spacing, they are ''relaxed'' towards the Boltzmann equilibrium distribution.These two steps applied to   which relates to the fluid density  are found to yield the Navier-Stokes equations for incompressible fluids (Chen and Doolen, 1998).The two steps applied to   are found to yield the convection-diffusion equation for energy density  which relates to the temperature  (Krüger et al., 2017).
The combined movement and collision steps can be written as and where   is the velocity vector pointing between adjacent lattice sites (i.e.   equals one lattice spacing orthogonally or diagonally in a square or cubic lattice), and    and    are the collision terms which redistribute number densities between directions due to collisions.These collision terms can be calculated by an efficient method of relaxing the distributions towards the equilibrium distribution due to Bhatnagar, Gross and Krook 1954 using where   and   are the relaxation times for the two number densities   and   , and    and    are respectively the corresponding equilibrium distributions which can be calculated as a Taylor expansion to second order of the Boltzmann distribution, namely and where  is the macroscopic density,  is the macroscopic fluid velocity,  is the energy density,   are weights which depend on the lattice (to be specified shortly) and   = ∕ √ 3 is the speed of sound in the lattice where  = ∕ is the lattice speed (i.e. is the lattice spacing of the cartesian lattice and  is the time step).In the following, we use non-dimensionalized simulations and hence, we set  =  = 1 which leads to equilibrium distributions given by and The relaxation times in the above equilibrium distributions relate to the kinematic viscosity and the thermal diffusivity as where   is the kinematic viscosity and where   =  is the thermal diffusivity.So far, the above equations are valid for 2D or 3D.However, in the following numerical study, we will perform simulations in 2D.We use the so called D2Q9 lattice (2 dimensions, 9 lattice velocities) in which number densities move on a square lattice along the axes or in the diagonal directions.We define the velocity vectors in the D2Q9 lattice   as and the weights vector   for this lattice is given by i.e. a weight of 4/9 for the zero velocity vector (i.e. the velocity vector for no movement), 1/9 for directions along the cartesian axes, and 1/36 along diagonals (note that ∑    = 1).It is well known that the above choice leads to the Navier-Stokes equations for incompressible fluid flow (Chen and Doolen, 1998) where the macroscopic properties can be obtained from the number densities   and   using and where the energy density  relates to temperature through  =  ∕2 where  = 2 is the number of dimensions,  is temperature and  is the gas constant which we set to unity for our non-dimensionalized simulation runs in this paper.

Boussinesq buoyancy term
In order to model convection, we need to add a buoyancy force and we adopt the usual Boussinesq approximation where density variations have a fixed part plus a perturbation that linearly depends on temperature given by We assume gravity and hence, under the Boussinesq approximation, we need to add a force term of where  is the acceleration due to gravity and  is the gravitational force due to gravity.The above gravitational force can be added to the LBM by adding a forcing term to the Lattice Boltzmann equations for   .Namely, Eq. (1) including the gravitational force becomes where the Boussinesq buoyancy forcing term   is given by where  = || =   = − .Eq. ( 19) can be verified by calculating the change in velocity during one time step due to the Boussinesq term using Eq. ( 14).

Extremely nonlinear ''plastic'' rheology
In order to add the effect of a plastic rheology, the total viscosity denoted   is normally specified as a function of the strain rate (Tosi et al., 2015) such as where and where  * is the effective viscosity at high stress,   is the yield stress, ε is the shear strain rate, and  is a power relating to the nonlinear rheology which is normally  = 1 for plasticity and  = (1 − 1∕) for dislocation creep.With the above formula, one could set  > 1 to model a more strongly nonlinear strain rate weakening material than plastic rheology with  = 1.In this work, we consider setting  to any positive value as a generalization of standard plasticity, and we refer to cases with  ≫ 1 as extreme plasticity.In the above formula, one can see that the plasticity effect will be negligible if   ≫   .In other words, the plastic effect is negligible when ε ≪   assuming that the effective viscosity at high stress denoted  * is small, i.e.  * <   .
In the LBM, we can calculate the shear strain rate using (Chen et al., 2014) where   is the shear strain rate tensor given by where   ( − ) is the relaxation time at the last time step and    is the non-equilibrium component given by In past work (Mora and Yuen, 2018b), we demonstrated the ability of the LBM to model temperature dependent viscosity such as Reynolds viscosity of form with powers up to  = 15.Here, we will set   =  (i.e.no temperature or depth dependence of viscosity) and focus on demonstrating the ability of the LBM to model extremely strong plasticity.

Numerical experiment
In this section, we will apply the LBM described above to a two layer model with an upper low viscosity layer, and demonstrate the LBM's ability to model strong plasticity which leads to the emergence of plate tectonic like behavior.

Setting of viscosity and thermal diffusivity
We first utilize the definition of the Rayleigh number given by where  =   is the kinematic viscosity,  =   is the thermal diffusivity,  is the gravitational acceleration,  is the coefficient of thermal expansion,  is the temperature difference between the top and bottom of the layer, and  is the layer thickness.Now, given that the Prandtl number   is defined by and using Eq. ( 27), we can set the viscosity and thermal diffusivity to where and hence, we can calculate the two relaxation times   and   required by the LBM using Eqs.( 9) and (10).

Two layer model initialization
We initialize a 2 layer model where the upper layer has a lower viscosity than the lower layer by a factor of 50, analogous to an assumption that the upper mantle has a 50× lower viscosity than the lower mantle.This is achieved by first calculating the viscosity and thermal diffusivity based on a Rayleigh number of  = 4.1667 × 10 6 and   = 1000, assuming a single layer.The dimensions of the computational lattice in  and  are set to  = 512,  = 256.Hence, from Eqs. ( 29) and ( 30) and setting  = 0.41667,  = 0.1,  = 0.001, and layer thickness1  =  − 3 = 253, we obtain the viscosity   and thermal diffusivity   =  to be We then set the thickness of the upper layer to be approximately 20% of the total model thickness  (i.e.similar to the upper mantle thickness relative to the total mantle thickness) so we have the upper layer thickness  1 = 50 ≈ 0.2 × ( − 3).Finally, we set the viscosity of the lower layer to   2 =   = 0.4024, and the viscosity of the upper layer to be 50× smaller   1 =   ∕50 = 0.00805.These viscosity values then allow the relaxation times for all parts of the model to be calculated using Eqs.( 9) and (10).

Setting of plasticity parameters
We set the power  in Eq. ( 22) in the upper and lower layers to be  1 = 4 and  2 = 1, so the plasticity will be much stronger in the upper layer (i.e.extreme plasticity in our terminology).The value of  * is set to  * = 0.00102 which is significantly less than the viscosity of the upper layer (i.e. about 8× less), and the value of the yield stress   is calculated as follows where  = 1, 2 specifies which layer.The above formula means that when the strain rate exceeds 80% of the maximal strain rate, the plasticity effect will become significant.A summary of the model parameters is given in Table 1.

Boundary conditions
We use periodic boundary conditions in the -direction, no-slip boundary conditions at the base of the lower layer, and a free slip boundary condition at the top of the upper layer.The no-slip boundary condition at the base of the lower layer is achieved through the standard LBM bounce-back approach.Namely, at the base of the lower layer at  =   we set the upgoing number densities as the bounced back downgoing number densities, namely  4 =  3 ,  7 =  8 and  6 =  5 .At the top of the upper layer at  =   , we developed free slip in  boundary conditions based on the approach of Zou and He (1997) in which we solve for ,   and downgoing number densities  3 ,  5 and  8 in terms of the known number densities  1 ,  2 ,  4 ,  6 and  7 and the known -component of velocity   = 0.This development uses Eqs. ( 13) and ( 14) and the bounce-back of non-equilibrium parts as in Zou and He (1997) The boundary conditions for  are the standard LBM bounce-back boundary conditions at the upper and lower boundaries.So at the top of the model at  =   we have  3 =  4 ,  5 =  6 and  8 =  7 , and at the base of the model at  =   we have  4 =  3 ,  6 =  5 and  7 =  8 .These boundary conditions amount to thermally insulating boundary conditions.In the runs in the next section, we will initialize a constant temperature  =  0 in the lattice, and set the temperature of the upper and lower boundaries to be   =  0 −  ∕2 +   and   =  0 +  ∕2 +   , where  = 0.001 as per our equations to calculate the layer viscosities in Eq. ( 29) and  is small zero mean Gaussian distributed random fluctuations with standard deviation of 2.5% of  .These random fluctuations are required to help to initiate convection.

Numerical stability
The thermal LBM is stable for a large range of viscosities, Rayleigh numbers and Prandtl numbers as shown in Mora and Yuen (2018a), and it is stable in the present work.However, it should be noted that this is not always so.Instabilities -if they occur -are generally due to either the fluid velocity speed || becoming too large relative to the   (i.e. the Mach number being too high) which leads to the equilibrium distribution ''blowing up'', and/or the non-dimensional viscosity being too small.
Regarding the choice of viscosity   , the LBM community generally thinks in terms of the relaxation time   .For a value of   that is too small, the LBM will become unstable.From Eq. ( 9), one can see that   has a minimum value of 0.5 when the viscosity is zero.Hence, one must choose a value of   >    > 0.5 to ensure stability of the LBM where the value of    depends on the problem being solved.For example, in the problem described above, the lowest possible viscosity being the effective viscosity at high stress denoted  * = 0.00102 which means that the lowest possible value for the relaxation time in the model is   * ≈ 0.503, and our simulation was stable with such low values for   .This implies that for our problem with an extremely plastic rheology,    for stability is within the range 0.5 <    ≲ 0.503.
Finally, if stability is an issue for a given problem using the BGK Single Relaxation Time (SRT) LBM method described above, the LBM's stability can potentially be improved by using the Multiple Relaxation Time (MRT) LBM method (d 'Humières et al., 2002;Aslan et al., 2014) in which one calculates the relaxation to equilibrium in moment space.

Getting the time scale and converting to physical units
In the above, we chose values for the various quantities in lattice units such that we ensured stability without thought to the time scale implied with our choices.However, the duration of each time step can be converted to seconds from lattice units as follows.
We start by re-defining  as the number of meters per lattice space unit, and  as the number of seconds per lattice time unit (one time step) = amount of physical time per time step.Hence, we have the following relationships between lattice and physical units where the subscript  denotes that a quantity is in physical units, and the subscript  denotes that a quantity is in lattice units.Hence, from the above relationships, and assuming we have chosen a value for  = spatial resolution of the LBM in physical units, we can derive the value of  number of seconds corresponding to one time step as Note that the above calculation assumes that the values of     =     so if this is not so, then before applying Eq. ( 35), one must convert   and   to the desired physical values of   and   while leaving the Rayleigh number the same by adjusting   .
To illustrate, we will now calculate the time scale for the runs in this paper assuming the physical quantities are as follows   = 4.5 × 10 3 kg m −3   = 9.8 m s −2 ,   = 3 × 10 −5 (1∕  (36) which are similar to the parameters for the Earth's mantle, where the value for   above is back calculated from the other values above assuming our Rayleigh number of   =   = 4.1667 × 10 6 .The above values lead to a physical Prandtl number of    =   ∕  = 6.617 × 10 24 which is very high and is often assumed to be infinity in geodynamical simulations.The LBM calculations are dynamic and cannot be done at infinite Prandtl number.However, if the value of   is large enough, a dynamic solution such as by using the LBM will asymptotically approach the solution at an infinite Prandtl number.
In the specifications of our numerical experiment above, we assume   = 10 3 although previous work suggests that the LBM can achieve higher Prandtl numbers of at least   = 10 4 (Mora and Yuen, 2018a).
Assuming the physical values for the earth given in Eq. ( 36), we can now rescale to obtain the physical time step  using Eq. ( 35) as follows.First, we modify the value of   ,   and   to obtain an equivalent set of values where   ′ is a new value for gravitational acceleration in lattice units re-scaled such that the Rayleigh number is identical but with lattice values for   ′ and   ′ which are identical to the desired physical values given in Eq. (36).
To obtain the time scale, the values for   ,   ′ ,   ,   and  are substituted into Eq.( 35).In our simulations, we have   = 2.9 × 10 6 m so for our lattice grid size of  = 256 we have   = ( − 3) = 253 ⇒  =   ∕  = 1.146×10 4 m.Applying Eq. ( 35 Hence, the  = 700 thousand time step simulations in the following corresponds to a total equivalent physical simulation time of   =  = 1.248 billion years.We note that by adjusting the values of   and   while maintaining the same Rayleigh number, we could achieve the same physical simulation time   with a greater or smaller time step  and hence, using a lower or higher number of time steps.

Results
We ran two simulations for the two layer model specified above.The first used only a viscous rheology, so we set   =   =  where   and   were as specified in Table 1.In the second, we applied Eq. ( 20) to calculate the total viscosity with the plasticity calculated using Eq. ( 22) and used the values for   ,   ,  * ,   and  specified in Table 1, so ''extreme plasticity'' in the upper layer.
Fig. 1 shows snapshots of the magnitude of the velocity || for the extreme plasticity case specified in the preceeding section (left column) and for the linear viscosity case where   =   (right column).This plot enables the general style of convection to be seen for each case.For the case of extreme plasticity (left column), one observes convection cells in the upper layer initially ( = 100 K).However, later in the simulation, the convection in the upper layer tends to consist of narrow horizontal channels along the upper edge of the model (see  = 300 K, 600 K and 700 K time steps), separated by upwelling or downwelling plumes.This can be considered as plate tectonic like behavior in that convection is mainly confined to narrow horizontal zones near the top of the upper layer and to vertical plumes.However, sometimes this plate tectonic like behavior is disrupted such as at  = 500 K time steps where one sees some convection in a broader zone in the upper layer and some evidence for convection cells within the upper layer.For the viscous case with no plasticity (right column), one tends to observe either convection cells within the upper layer, broad horizontal zones and plumes within the upper layer, or a combination of these two modes of convection.Note that we obtain plate like behavior with downwelling plumes rather than plate tectonics with subduction zones because we are focussed on modeling extreme plasticity whereas plate tectonics with subduction zones requires stiff plates which can be achieved by the addition of temperature dependent viscosity (Tackley, 2000).
Fig. 2 shows snapshots of vector plots of  * = (  , −  ) which show the flow pattern.One can clearly see on the plots for the case of extreme plasticity (left column) that there is plate like flow in the upper layer at  = 300 K, 600 K and 700 K time steps, with most of the flow confined to a narrow layer at the top of the model, or to upwelling plumes (c.f.spreading centers) and downwelling plumes (c.f.subduction zones).Again, one observes a different pattern of convection for the viscous case (right column) with convection in the upper layer being more distributed within the layer as smaller convection cells, or horizontally flowing broad regions within the upper layer.
Fig. 3 shows snapshots of temperature  for the extreme plasticity case (left column) and linear viscosity case (right column).For the extreme plasticity case, one can see the detailed form of the upwelling plumes below the spreading centers and downwelling plumes below the convergence zones at  = 300 K, 600 K and 700 K. Again, snapshots for the linear viscosity case show convection tends to occur throughout the upper layer.
Regarding computational cost of the 700 K time step = 1.25 Gy simulations above, the extremely plastic simulation took 17.4 min on 16 CPU 2.2 GHz Xeon cores and the viscous simulation took 13.8 min.The LBM kernel calculations ran at 3.45 GFlops = 78% of peak performance per core, and the total LBM algorithm including streaming and inter-core communications ran at 1.93 GFlops = 44% of peak performance per core.The computer code runs in 2D or 3D and was written in Python and parallelized using MPI with the LBM kernel calculated using optimized C code, and has been shown to have linear weak scaling to thousands of cores (Mora et al., 2020(Mora et al., , 2023)).Note that further optimizations to the streaming and communications steps may be possible which would further improve the performance of the total LBM algorithm.

Use of geodynamic simulations in AI
Machine learning applied to geodynamics data is still in its infancy, however applications have been put forward to investigate potential practical applications.At the level of whole mantle convection, thanks to mantle tomography, present Earth's detailed velocity and density structure is relatively well known, however the past history and detailed interiors of terrestrial planets, either the ones in our solar system or newly discovered exo-planets, can only be indirectly inferred.Due to the complexities of the differential equations driving mantle convection, coverage of the entire parameter space of mantle convection is infeasible, even in 2D.Machine learning has been shown to be able to at least partially overcome this computational bottleneck by running end-member cases of large scale 2D or 3D simulations and then extracting potential observable parameters (e.g.1D seismic profile for terrestrial planets, thermal emission and atmosphere composition for exo-planets).Examples of this approach are in Agarwal et al. 2020.
Similarly, machine learning has been used to process large 3D models and extract snapshots of convection models as training, testing and verification samples.For example, Shahnas et al. (2018) used Support Vector Machines (SVMs) to predict mantle density anomalies and mantle flow patterns, illustrating how machine learning can be used to add time dependency to geodynamic modeling.Snapshots of convection models can also be used to analyze the effect of varying mineral physics data, which can rapidly expand the parameter space of the models to be explored.
Supervised learning techniques that can be used to analyze full convection models vary from Convolutional Neural Networks (CNNs) (Morra et al., 2020), which are already widely used to analyze images in geosciences, as well as Support Vector Machines.For example, recently efforts have been made to extract fundamental parameters such as Prandtl and Rayleigh numbers from a CNN applied to LBM models of convection (Boroumand et al., 2023) where 200 snapshots of 88 convection simulations were used to train and test the Neural Network based classifier.In a future not too remote, we can see how Generative Artificial Intelligence will be used to predict solutions of mantle dynamics models, as recent attempts in atmospheric convection (Mooers et al., 2020) and convection in general (Jiang and Farimani, 2020) have shown to be possible.
While all these approaches are promising and give us a snapshot on the future of geodynamics modeling, the present challenge is represented by the large number of samples of training set and by the high level of quality, detail and reliability of the training set.In this sense, LBM represents a unique opportunity for geodynamics research to produce the necessary datasets, due to the capability to almost perfectly scale to large number of processors (Mora et al., 2020) with very high performance (Mora et al., 2023), and to model at the same time high as well as low Prandtl numbers (Mora and Yuen, 2018a).Low Prandtl number models are necessary to model the formation of terrestrial planets, when all or part of the mantle was molten, and to investigate the transition to the present solid state convection, which is the engine of plate tectonics.

The plate tectonic-like behavior defined as complexity theory
Another application of large scale 2D and 3D models is to finally elucidate some of the mysteries of plate tectonics, which require multiple space and temporal scales to be accurately reproduced.Major examples in geodynamics are the feedback at different scales between the megathrust scale and the subduction zone scale, which has seen major recent strides in numerical modeling (Van Dinther et al., 2013), as for example in Sobolev and Muldashev (Sobolev and Muldashev, 2017), which reflects advances in laboratory models, Corbi et al. (2020), Johnson et al. (2021).
Also the long term geodynamic scale requires high resolution models in 3D, until now not available.For example Fractal Plate tectonics (Sornette and Pisarenko, 2003) is likely an expression the complexity of the interplay between the large scale mantle flow and the brittle rupture of the surface (Bird, 2003) as suggested by time dependent analysis of the last 200 Myrs of plate boundary evolution (Morra et al., 2013).Large scale models have shown that part of the scaling between mantle flow and lithospheric deformation reflects the mantle flow patterns (Mallard et al., 2016) but the fractal scaling across the scales is still a mystery (Bercovici et al., 2015).

Computational performance
We have not directly compared the computational performance of the LBM to that of other approaches which would require comparisons of specific simulations within the range of applicability of both methods.For example, many classical approaches are developed assuming infinite Prandtl number while the LBM can model a wide range of Prandtl numbers such as 1 − 10 4 as well as Rayleigh numbers up to over  ∼ 10 12 .Recent work Mora et al. (2023) indicates that the LBM is able to achieve extremely high performance on up to thousands of cores with the computational kernel running at around 4-5 GFlops/core, and the overall algorithm including the streaming step and MPI communications running at 1.3-3 Gflops per core.The exact computational performance within these ranges depending on number of cores being used and hardware details such as cache sizes combined with memory and infiniband bandwidth of the HPC cluster.Whether the LBM or classical methods are more efficient for any specific simulation where both are applicable is yet to be studied.However, given the example in this paper where we modeled a highly nonlinear rheology case for 1.25 Gy in 17.4 min on just 16 cores, we believe that the LBM performs extremely well for many problems in geodynamics.
Furthermore, given the LBM's ability to model highly nonlinear rheologies, sharp viscosity contrasts, and a wide range of Rayleigh & Prandtl numbers, we believe the that at least in some cases, the LBM goes beyond the capabilities of classical methods.Hence, we believe that in many problems -where both methods can be applied -the LBM will be computationally efficient and/or superior to classical methods.Furthermore, the LBM also offers extended geodynamical modeling capabilities such as for simulations at extremely high Rayleigh numbers and low Prandtl numbers and efficiently modeling nonlinear rheologies.We note however that at high Prandtl numbers where the dynamics asymptotically approaches that of an infinite Prandtl number, the LBM may require a large number of time steps to initiate convection.As such, for simple rheologies (e.g.viscous) at extremely high Prandtl numbers whose dynamics asymptotically approached that of infinite Prandtl number, classical methods may be more efficient than the LBM although this is yet to be studied.
Unfortunately, to our knowledge, although a direct comparison of the performance of several codes have been done for the geodynamo (Hiroaki et al., 2016), no such analysis has been performed for generally used mantle convection codes, such as Aspect (Bangerth et al., 2020), Underworld (Mansour et al., 2020), or GAIA (Hüttig and Stemmer, 2008).For these reasons we can only mention that while the first two of these codes use a sparse matrix-based solver, and GAIA a matrix multiplication one, our approach is matrix-free and therefore promises to ultimately be able to scale more successfully for large number of processors.Furthermore, our code works in both 2D and 3D.Existing scaling tests exist for some of the above mentioned software, e.g. for Aspect, scaling performed for the XSEDE cluster (Anon,0000) show that linear scaling is lost at about 100 cores for 2.7 × 10 5 Degrees of Freedom (DOFs) and at 1000 cores for about 10 7 DOFs, so a model resolution of order a few hundred cubed.Similar results are have been found for Underworld, where linear scaling is found to around 1000 cores for a 3D mesh of 256 3 so 16.8 million model size, estimated equivalent to 7×10 7 DOFs (Giordani, 2021).While a direct performance comparison would require running exactly the same problem, with the same resolution and for the same scaled geodynamical time, our LBM implementation displays weak linear scaling to 2048 cores (the largest number available to us) for up to 64 million computing points where further details and performance analysis will be given in a publication in preparation (Mora et al., 2023).

Conclusions
The Lattice Boltzmann Method (LBM) is a semi-microscopic method in which particle number densities move and collide on a discrete lattice which has been shown to model the incompressible Navier-Stokes equations in the macroscopic limit.We have presented the thermal LBM as an alternate means to classical numerical solutions to PDE's for geodynamic simulations of mantle convection in the Earth or exoplanets.The advantage of the LBM is that it is able to handle large viscosity ranges and contrasts, temperature dependent rheologies and highly nonlinear rheologies such as ''extreme plasticity'', which endows mantle convection with a connection to the past.We demonstrate the method by modeling a two layer system with an upper low viscosity layer with a highly nonlinear rheology with    ∝ ε− where  = 4.We find that, as expected, introduction of this highly nonlinear rheology leads to the emergence of plate tectonic like behavior where emergence is as defined by complexity theory (Nicolis and Prigogine, 1989).We therefore propose the LBM provides a valuable alternative to classical numerical solutions to PDE's for modeling the dynamics of  the earth and exoplanets.The LBM's ability to model a wide range of viscosities, multiphase flow, strong viscosity contrasts, high Rayleigh numbers, a wide range of Prandtl numbers and highly nonlinear rheologies such as ''extreme plasticity'' with    ∝ ε−4 as shown here, potentially allows it to model the dynamics of the earth and exoplanets from the early lava world stage through to plate tectonics or other regimes.Our code is available to researchers upon request for collaborative research and in time will become open source, and moreover, it is parallelized using MPI with linear weak scaling, and hence, has the potential to run on 10 4 to 10 5 cores, which can make use of today's era of exascale computing.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Snapshots of the magnitude of velocity || at various time steps with extreme plasticity (left) and for the viscous case -i.e.no plasticity (right).Note that the values are clipped at half of ||  so dark red shows regions where || > 0.5||  .
P.Mora et al.

Fig. 2 .
Fig. 2. Snapshots showing a vector plot of velocity  * = (  , −  ) at various time steps with extreme plasticity (left) and for the viscous case -i.e.no plasticity (right).Note that the color scale shows either   if |  | > |  | so red = right and blue = left, or −  otherwise so red = up and blue = down.

Fig. 3 .
Fig. 3. Snapshots of temperature  at various time steps with extreme plasticity (left) and for the viscous case -i.e.no plasticity (right).Cold downwellings can be seen as blue colors, and upwelling plumes can be seen in red colors.

Table 1
Model parameters.