We present 2D simulations using the Lattice Boltzmann Method (LBM) of a fluid in a rectangular box being heated from below, and cooled from above, with a Rayleigh of$Ra\mathrm{<mo>=</mo>5}\times {10}^7$ , similar to current estimates of the Earth's mantle, and a Prandtl number of 1,000 which assures simulations lie within the non-inertial regime relevant to geodynamics. Simulations are conducted with two different temperature dependent viscosities, Reynolds-type temperature dependence with$\nu \mathrm{<mo\; stretchy="false">(</mo>}T\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\nu }_0\exp \mathrm{<mo\; stretchy="false">(</mo>}T\mathrm{<mo>/</mo>}T\ell \mathrm{<mo\; stretchy="false">)</mo>}$ , and according to the Arrhenius law with$\nu \mathrm{<mo\; stretchy="false">(</mo>}T\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\nu }_0\exp \mathrm{<mo\; stretchy="false">(</mo>}E\mathrm{<mo>/</mo><mo\; stretchy="false">(</mo>}RT\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo>=</mo>}{\nu }_0\exp [b(T_0\mathrm{<mo>/</mo>}T-\mathrm{1<mo\; stretchy="false">)</mo>}/\mathrm{<mo\; stretchy="false">(</mo>}{T}_0\mathrm{<mo>/</mo>}{T}_u-\mathrm{1<mo\; stretchy="false">)</mo>}]$ . We conduct numerical experiments for exponent b  in the range 7 through 15. For both cases, we achieve stagnant-lid type convection. For the case of Reynolds temperature dependent viscosity, we obtain increasingly pulse-like plumes as $b$  is increased where plume-head bubbles of hot upwelling material push upwards through pre-weakened pathways to the surface of the model. In contrast, the Arrhenius temperature dependent viscosity leads to long narrow plumes separated by long downwelling chutes that push through the entire model. These results demonstrate that the LBM can model different modes of convection for different temperature dependence of viscosity, and point to the need to accurately define the temperature dependence of viscosity in convection problems such as geodynamics.

Keywords: Lattice Boltzmann Method, natural convection, Reynolds and Arrhenius temperature dependent viscosity, Boussinesq approximation, Rayleigh number, Prandtl number

Simulation and visualization of geodynamical processes is a scientific challenge due to spatial and time scales spanning many orders of magnitude and highly nonlinear phenomena.^1-3 Research to date has focused largely on using direct methods to solve the partial differential equations describing the conservation of mass, momentum, and energy in the macroscopic domain.^4-7 An alternative approach has arisen in the last few decades in which one simulates a discrete Boltzmann equation on a lattice that describes the movement and collision of distributions of particles, rather than solving the macroscopic conservation equations themselves. These Lattice Boltzmann Method have been proven to yield the Navier-Stokes equations in the macroscopic limit,⁸ by using a Chapman-Enskog expansion of the number densities in terms of macroscopic velocity, just as was done for the continuous Boltzmann equation.⁹

Lattice Boltzmann methods have their origins in Lattice Gas Automata which were first proven by Frisch et al.,¹⁰ to yield the Navier-Stokes equations in the macroscopic limit. These initial models were unconditionally stable and conserved mass, momentum and energy perfectly. However, they were computationally expensive with averaging needed over space to obtain the macroscopic equations and furthermore, costly calculations were required to compute the collision term. However, since the development of an efficient method via relaxation to calculate the collision.^8,11,12 Research and applications of the Lattice Boltzmann method have recently undergone an explosion.¹³ Numerous studies have been conducted of thermal convection such as He at al.,¹⁴ Guo et al.,¹⁵ Shan,¹⁶ Wang et al.,¹⁷ Arun et al.¹⁸ and Huber¹⁹ models melting and convection applied to geological problems. Many papers focus on model validation and study Rayleigh-Bénard convection, with some investigating the behavior of plume growth in high Rayleigh number flows up to$Ra\mathrm{<mo>=</mo>4.9}\times {10}^9$ .²⁰ Given the current estimates of the Rayleigh number of the Earth's mantle of order ${10}^7-{10}^8$ ,¹ it seems that the LBM offers an alternate possibility to simulate and study mantle convection.²¹ Outstanding scientific advances have been made on understanding the style and platform of mantle convection through numerical simulations via macroscopic partial-differential equations.²²⁻²⁵ Up to now, Lattice Boltzmann Methods have not been widely used in the mantle convection community and may offer an alternate approach to studies of mantle convection with more complex rheology. The strengths of the LBM is that this method is relatively simple to code and is also easy to parallelize and put on accelerators such as GPU's, it is capable of modeling multi-phase fluid flow and phase transitions such as melting or polymorphic solid-state transitions, and it does not require interfacial tracking to model multiphase flow. Hence, surface tension and other interface effects are automatically handled by the method.²⁶ Huang et al.,¹³ provide an extensive overview of multiphase Lattice Boltzmann methods which are currently being researched and applied.

In this paper we present simulations of thermally driven convection with strong temperature dependent viscosity with the aim of demonstrating the potential of the Lattice Boltzmann Method to be applied to problems in geodynamics, and studying the effect of strongly temperature dependent viscosity on plume dynamics.

The Lattice Boltzmann Method (LBM) involves simulating distributions of particles moving and colliding on a discrete lattice. In one time-step, the particles can move by one lattice spacing along the orthogonal axes, or along diagonals. We use the standard notation in the LBM denoted DnQm for a simulation in $D\mathrm{<mo>=</mo>}n$ dimensions, and with $Q\mathrm{<mo>=</mo>}m$ velocities on the discrete lattice. In the following, we restrict ourselves to 2D and use the $D2Q9$  Lattice Boltzmann lattice arrangement shown in Figure 1. Appendix A shows the nomenclature used in this paper. In the lattice, we define $f_{\alpha }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$  as the number density of particles moving in the $\alpha$ -direction where the $Q\mathrm{<mo>=</mo>9}$  velocities are given by$e_{\alpha }\mathrm{<mo>=</mo><mo\; stretchy="false">[</mo><mo\; stretchy="false">(</mo>0,0<mo\; stretchy="false">)</mo>,<mo\; stretchy="false">(</mo>1,0<mo\; stretchy="false">)</mo>,<mo\; stretchy="false">(</mo>0,1<mo\; stretchy="false">)</mo>,<mo\; stretchy="false">(</mo>}-\mathrm{1,0<mo\; stretchy="false">)</mo>,<mo\; stretchy="false">(</mo>0,}-\mathrm{1<mo\; stretchy="false">)</mo>,<mo\; stretchy="false">(</mo>1,1<mo\; stretchy="false">)</mo>,<mo\; stretchy="false">(</mo>}-\mathrm{1,1<mo\; stretchy="false">)</mo>,<mo\; stretchy="false">(</mo>}-\mathrm{1,}-\mathrm{1<mo\; stretchy="false">)</mo>,<mo\; stretchy="false">(</mo>}-{\mathrm{1,1<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>}}^T$ . This choice means that $e_0$  is the zero velocity vector and so represents stationary particles, and $e_{\alpha }\mathrm{<mo>=</mo>}-e_{\alpha +2}$ for$\alpha \mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>1,2,5,6<mo\; stretchy="false">)</mo>}$ . The lattice is unitary so the lattice spacing and time step are$\Delta x\mathrm{<mo>=</mo>}\Delta t\mathrm{<mo>=</mo>1}$ .

The thermal Lattice Boltzmann Method involves two steps, movement and collision, of two distribution functions.¹⁴ If we wish to model thermal-chemical convection, we just add in another distribution function, representing the chemical component. In this paper on thermal convection, we therefore define two distribution functions, $f_{\alpha }$ representing the mass density of particles moving in the $\alpha$ -direction, and $g_{\alpha }$  representing the energy density of particles moving in the $\alpha$ -direction. The evolution equations encompassing the two steps are given by
$f_{\alpha }\mathrm{<mo\; stretchy="false">(</mo>}x+e_{\alpha }\Delta t\mathrm{,}t+\Delta t\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}{f}_{\alpha }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}+\Delta f_{\alpha }^c\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>},$                             (1)
and
$g_{\alpha }\mathrm{<mo\; stretchy="false">(</mo>}x+e_{\alpha }\Delta t\mathrm{,}t+\Delta t\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}{g}_{\alpha }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}+\Delta g_{\alpha }^c\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>},$                            (2)

where $\Delta f_{\alpha }^c\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$  is the collision term and represents the redistribution of particle number densities at lattice site $\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$ due to collisions. The collision term can be calculated exactly, but this is computationally expensive and as such, is rarely done. Alternatively, the collision term can be calculated by the BGK method in which case, the distributions relax to the equilibrium distribution.^8,12 The BGK method is computationally efficient and gives satisfactory results provided the distributions are not too far from equilibrium. The BGK collision term is given by
$\Delta f_{\alpha }^c\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}\left(\frac{1}{{\tau }_f}\right)\mathrm{<mo\; stretchy="false">(</mo>}{f}_{\alpha }^{eq}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}-f_{\alpha }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>},$                                    (3)
and
$\Delta g_{\alpha }^c\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}\left(\frac{1}{{\tau }_g}\right)\mathrm{<mo\; stretchy="false">(</mo>}{g}_{\alpha }^{eq}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}-g_{\alpha }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>},$                                    (4)

where $f_{\alpha }^{eq}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$  is used to denote the equilibrium distribution of $f_{\alpha }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$  and is given by
$f_{\alpha }^{eq}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}\rho w_{\alpha }\left[1+C_1\frac{\mathrm{<mo\; stretchy="false">(</mo>}{e}_{\alpha }\cdot u\mathrm{<mo\; stretchy="false">)</mo>}}{c^2}+C_2\frac{{\mathrm{<mo\; stretchy="false">(</mo>}{e}_{\alpha }\cdot u\mathrm{<mo\; stretchy="false">)</mo>}}^2}{c^4}+C_3\frac{{\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}}^2}{c^2}\right],$  (5)
and $g_{\alpha }^{eq}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$  is used to denote the equilibrium distribution of $g_{\alpha }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$  and is given by
$g_{\alpha }^{eq}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}\rho \epsilon w_{\alpha }\left[1+C_1\frac{\mathrm{<mo\; stretchy="false">(</mo>}{e}_{\alpha }\cdot u\mathrm{<mo\; stretchy="false">)</mo>}}{c^2}+C_2\frac{{\mathrm{<mo\; stretchy="false">(</mo>}{e}_{\alpha }\cdot u\mathrm{<mo\; stretchy="false">)</mo>}}^2}{c^4}+C_3\frac{{\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo\; stretchy="false">)</mo>}}^2}{c^2}\right],$  (6)

where $c\mathrm{<mo>=</mo>}\Delta x\mathrm{<mo>/</mo>}\Delta t\mathrm{<mo>=</mo>1}$  is the lattice speed and the constants $C_1$  through $C_3$  are given by $C_1\mathrm{<mo>=</mo>3}$ , $C_2\mathrm{<mo>=</mo>9<mo>/</mo>2}$  and $C_3\mathrm{<mo>=</mo>}-\mathrm{3<mo>/</mo>2}$ . The speed of sound in the $D2Q9$  LBM lattice is given by $c_s\mathrm{<mo>=</mo>1<mo>/</mo>}\sqrt{3}$ .

In Equation (3), the value of ${\tau }_f$  is a relaxation time which relates to the kinematic viscosity ${\nu }_f$ through
${\tau }_f\mathrm{<mo>=</mo>}{\nu }_f\mathrm{<mo>/</mo><mo\; stretchy="false">(</mo>}{c}_s^2\Delta t\mathrm{<mo\; stretchy="false">)</mo>}+0.5.$           (7)
Similarly, ${\tau }_g$ in Equation (4) is a relaxation time which relates to the thermal viscosity ${\nu }_g$  (and hence, the thermal diffusivity$\kappa \mathrm{<mo>=</mo>}{\nu }_g$ ) through
${\tau }_g\mathrm{<mo>=</mo>}{\nu }_g\mathrm{<mo>/</mo><mo\; stretchy="false">(</mo>}{c}_s^2\Delta t\mathrm{<mo\; stretchy="false">)</mo>}+0.5\mathrm{<mo>=</mo>}\kappa \mathrm{<mo>/</mo><mo\; stretchy="false">(</mo>}{c}_s^2\Delta t\mathrm{<mo\; stretchy="false">)</mo>}+0.5.$                                (8)

Figure 1The $D2Q9$ lattice.

In Equations (5) and (6) for the equilibrium distributions, the weighting scalars $w_{\alpha }$  are given by $w_0\mathrm{<mo>=</mo>4<mo>/</mo>9}$  for $\alpha \mathrm{<mo>=</mo>0}$  (stationary particles), $w_{\alpha }\mathrm{<mo>=</mo>1<mo>/</mo>9}$  for $\alpha \mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>1,2,3,4<mo\; stretchy="false">)</mo>}$ (particles travelling along the two cartesian axes), and $w_{\alpha }\mathrm{<mo>=</mo>1<mo>/</mo>36}$  for $\alpha \mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>5,6,7,8<mo\; stretchy="false">)</mo>}$  which are the particles travelling diagonally.

The macroscopic properties, density$\rho$ , velocity $u$ and temperature $T$ relate to the distribution functions $f_{\alpha }$  and $g_{\alpha }$  through
$\rho \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}{\displaystyle \sum _{\alpha }}{f}_{\alpha }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{,}$                                                             (9)
$P\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}\rho u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}{\displaystyle \sum _{\alpha }}{f}_{\alpha }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}{e}_{\alpha },$                                            (10)
and
$\rho \epsilon \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}{\displaystyle \sum _{\alpha }}{g}_{\alpha }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>},$                                                             (11)

where $\rho \mathrm{<mo>=</mo>}\rho \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$ is the macroscopic density, $P\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$ is the momentum density, and $\epsilon \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$  is the internal energy density which relates to temperature through
$\rho \epsilon \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}\rho DRT\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>2}\Rightarrow T\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}2DR\epsilon \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>},$               (12)

where $R$  is the gas constant and $D$  is the number of dimensions. In the following, we restrict ourselves to two dimensions ($D\mathrm{<mo>=</mo>2}$ ) and we use units such that$R\mathrm{<mo>=</mo>1}$ . Hence, we have$T\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\epsilon \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$ .

Buoyancy term - the Boussinesq approximation
In the following, we assume the Boussinesq approximation for the buoyancy term. Namely, density variations have a fixed part, plus a term that is linearly related to temperature as follows
$\rho \mathrm{<mo>=</mo>}{\rho }_0\mathrm{<mo\; stretchy="false">(</mo>1}-\beta \Delta T\mathrm{<mo\; stretchy="false">)</mo>},$                                                                  (13)
where $\beta$  is the coefficient of thermal expansion. Hence, if the body force is the gravitational force given by$G\mathrm{<mo>=</mo>}\rho g$ , we need to add a body force term to the Boltzmann equation given by equation (1) that is equal to
$G\mathrm{<mo>=</mo>}-\rho \beta \Delta Tg.$                                                                        (14)
This can be achieved by adding another term to equation (1). Namely,
$f_{\alpha }\mathrm{<mo\; stretchy="false">(</mo>}x+e_{\alpha }\Delta t\mathrm{,}t+\Delta t\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}{f}_{\alpha }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}+\Delta f_{\alpha }^c\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}+\Delta f_{\alpha }^b\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>},$ (15)
where the term $\Delta f_{\alpha }^b\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$  is the Boussinesq buoyancy forcing term which can be computed by setting $u\mathrm{<mo>=</mo>}G\Delta t$ into the equilibrium distribution given by Equation (5). Hence, to first order we have
$\Delta f_{\alpha }^b\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}\left(\frac{1}{{\tau }_b}\right)\mathrm{<mo\; stretchy="false">(</mo>}{w}_{\alpha }\mathrm{3<mo\; stretchy="false">(</mo>}{e}_{\alpha }{\mathrm{<mo\; stretchy="false">)</mo>}}_{z}G\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>},$                                      (16)
where $G\mathrm{<mo>=</mo><mo>|</mo>}G\mathrm{<mo>|</mo>}$ is the magnitude of the vertical gravitational force given by Equation (14), and ${\tau }_b$  is a relaxation time-like term that ensures the thermal buoyancy term has the same units as the other forcing terms of density. In this paper, we set ${\nu }_b\mathrm{<mo>=</mo>0.1}$  so ${\tau }_b\mathrm{<mo>=</mo>0.8}$ (cf. Equation (7)).

In the following, we conduct experiments in a rectangular domain of size $L_x\times L_z$  with a cold upper lid with temperature $T_u$  and a hot lower base with temperature$T_{\ell }$ . The Rayleigh number is a dimensionless number whose magnitude relates to the vigour of convection and is defined as follows
$Ra\mathrm{<mo>=</mo>}Gr\cdot Pr\mathrm{<mo>=</mo>}\frac{g\beta \mathrm{<mo\; stretchy="false">(</mo>}{T}_{\ell }-T_u\mathrm{<mo\; stretchy="false">)</mo>}{L}^3}{\kappa \nu },$                                             (17)
where Gr is the Grashof number relating buoyancy to viscous forces and Pr is the Prandtl number relating momentum and thermal diffusivity defined as follows
$Gr\mathrm{<mo>=</mo>}\frac{g\beta \mathrm{<mo\; stretchy="false">(</mo>}{T}_{\ell }-T_u\mathrm{<mo\; stretchy="false">)</mo>}{L}^3}{{\nu }^2},$                                                             (18)
$Pr\mathrm{<mo>=</mo>}\frac{\nu }{\kappa },$                                                                                                (19)
where $g$  is the acceleration due to gravity, $T_{\ell }$  is the temperature of the lower heated base, $T_u$  is the temperature of the upper cold lid, $v$  is the kinematic viscosity, $\beta$  is the coefficient of thermal expansion, $\kappa$  is the thermal diffusivity, and $L\mathrm{<mo>=</mo>}{L}_z$  is the vertical length scale between the lower (hot) and upper (cold) plates. Therefore, in the LBM, we can define the Rayleigh number as
$Ra\mathrm{<mo>=</mo>}\frac{g\beta \Delta T{L}^3}{{\nu }_g{\nu }_f},$                                                                        (20)
and
$Pr\mathrm{<mo>=</mo>}\frac{{\nu }_f}{{\nu }_g}$                                                                  (21)
where $g$ the acceleration due to gravity is, $\Delta T$  is the temperature range, $\beta$  is the coefficient of thermal expansion, ${\nu }_g$  is the thermal viscosity so $\kappa \mathrm{<mo>=</mo>}{\nu }_g$  is the thermal diffusivity, and $L\mathrm{<mo>=</mo>}{L}_z$  is the length scale. In the simulation shown in this paper, we have$L_x\times L_z\mathrm{<mo>=</mo>1024}\times 256$ . Hence,$L\mathrm{<mo>=</mo>}{L}_z\mathrm{<mo>=</mo>255}$ ,$g\mathrm{<mo>=</mo>1}$ ,$\beta \mathrm{<mo>=</mo>0.1},$ $\Delta T\mathrm{<mo>=</mo>}{T}_u-T_{\ell }\mathrm{<mo>=</mo>0.001}$  and with ${\nu }_f$  and ${\nu }_g$  set according to the desired Rayleigh number and Prandtl number of the flow. In the following, we conduct simulations with a Rayleigh number of $Ra\mathrm{<mo>=</mo>5}\times {10}^7$  and with a high Prandtl numbers of$Pr\mathrm{<mo>=</mo>1,000}$ .  The Rayleigh number relates to the vigour of buoyancy driven flow known as natural convection. Below a critical value$R{a}_{critical}\mathrm{<mo>=</mo>1708}$ , heat transfer is mainly by conduction, and above the critical value, heat transfer is mainly by convection. The value of the Rayleigh number of $5\times {10}^7$  was chosen to be in the mid-range of current estimates of the Rayleigh number of the Earth’s mantle which are of order ${10}^7–{10}^8{.}^1$ At this high Rayleigh number, we are in the regime of vigorous natural convection. At the high Rayleigh number of $5\times {10}^7$ and a Prandtl number of 1000, we are in the non-inertial regime where a series of narrow plumes are expected to form with mushroom heads and narrow necks, separated by narrow downwelling chutes.²⁷

It is well known that the viscosity of fluids depends on the temperature of the fluid, with viscosity decreasing with increasing temperature. Two commonly used models of temperature dependent viscosity are the exponential model first proposed by Reynolds and equivalent to the first order Frank-Kamenetsky approximation²⁸ and the Arrhenius model. The Reynolds' model is an empirical model that usually works for a limited range of temperatures. In this model, the viscosity depends on temperature according to the following formula
$\nu \mathrm{<mo\; stretchy="false">(</mo>}T\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\nu }_0\exp \mathrm{<mo\; stretchy="false">(</mo>}-bT\mathrm{<mo\; stretchy="false">)</mo>}.$                                                              (22)
The Arrhenius model assumes that the fluid flow obeys the Arrhenius equation for molecular kinetics and has the form
$\nu \mathrm{<mo\; stretchy="false">(</mo>}T\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\nu }_0\exp \left(\frac{E}{RT}\right),$                                                             (23)
where $R$  is the universal gas constant and $E$  is the "activation energy". In the following section, we do calculations with each model, and compare how these different models affect plume dynamics.

The Reynolds temperature dependence of viscosity is an empirical relationship that has a limited range of validity. It can be considered as a linearized approximation of the Arrhenius model. In the following, we test for various coefficients $b$  to determine when significant variations occur to the mode of convection.  In these experiments, we will set the the temperature of the upper lid to $T_u\mathrm{<mo>=</mo>}-0.0005$ and that of the lower base to $T_l\mathrm{<mo>=</mo>0.0005}$ so $T_0\mathrm{<mo>=</mo>0}$ . We set $b$  values such that the viscosity in the bulk of the model is lower by a factor of $\exp \mathrm{<mo\; stretchy="false">(</mo>}-b\mathrm{<mo\; stretchy="false">)</mo>}$  relative to the viscosity at the upper lid. Namely, we model $\nu \mathrm{<mo\; stretchy="false">(</mo>}T\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\nu }_0\exp \mathrm{<mo\; stretchy="false">(</mo>}-bT\mathrm{<mo>/</mo>}T\ell \mathrm{<mo\; stretchy="false">)</mo>}$ . Figure 2 shows snapshots at time step of 1,000,000, which is after the model achieves a steady state of convection for various coefficients $b$  ranging from $b\mathrm{<mo>=</mo>7}$ through$b\mathrm{<mo>=</mo>15}$ . These snapshots show that for low $b$ -values, the convection consists of narrow plumes separated by narrow downwellings progressing all the way through the layer. As the $b$ -value is increased, there is an increasingly pulse-like character of the upwelling plumes. This is due to the viscosity increase with temperature leading to bubbles of low viscosity light material pushing upwards through the model.  These bubbles tend to follow pre-weakened pathways of the preceding bubble trails. The upper cold thermal boundary layer is underlain by a hotter layer. However, the plots of the temperature profiles shown in Figure 3(A−D) show that this effect is minor for most values of the exponent$b$ . However, as the $b$ -value approaches $15$ , the mean temperature profile increases significantly with decreasing depth and is hottest just below the cold lid.

Figure 2 Snapshots of temperature in simulations using Reynolds temperature dependence of viscosity for cases of $b\mathrm{<mo>=</mo>7,}b\mathrm{<mo>=</mo>10,}b\mathrm{<mo>=</mo>12}$ and $b\mathrm{<mo>=</mo>15}$  (top – bottom).

A                                                          B

C                                                          D

Figure 3(A−D) Snapshots of mean temperature profiles with depth for simulations using Reynolds temperature dependence of viscosity for cases of $b\mathrm{<mo>=</mo>7,}b\mathrm{<mo>=</mo>10,}b\mathrm{<mo>=</mo>12}$ and $b\mathrm{<mo>=</mo>15}$ (top-left – bottom-right).

3D simulation results from the literature²⁹ show that linear upwellings rise from the Core-Mantle Boundary (CMB), and these quickly become quasi-cylindrical plumes which propagate upwards to the model surface. There is no pulsation observed in these plumes. This result is similar to the results in this paper obtained with the LBM when using weak temperature dependent viscosity, but contrast with the pulse-like behavior of plumes calculated by the LBM when temperature dependence is strong. The results of this paper suggest it is important to determine the strength of the temperature dependence of viscosity to model the Earth, as the form of plumes with strong temperature dependence is fundamentally different to the form of plumes when the temperature dependence is weak.

 In these runs we set the temperature of the upper lid to $T_u\mathrm{<mo>=</mo>0.0001}$  and the temperature of the lower base to $T_l\mathrm{<mo>=</mo>0.0011}$ . Figure 4 shows snapshots at a time step of 1,000,000 for the case of the Arrhenius model with different values of the activation energy, which is designed to be comparable the different $b$ -values of the Reynold's model for the upper Thermal Boundary Layer (TBL). Namely, for the upper TBL, we set $E$ such that the drop in viscosity in the bulk of the model compared to the viscosity at the upper lid is the same as that for $b$ -values of 7, 10, 12 and 15 in the Reynolds model.  This is achieved by setting the viscosity as
$\nu \mathrm{<mo\; stretchy="false">(</mo>}T\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}{\nu }_0\exp \mathrm{<mo\; stretchy="false">(</mo>}E\mathrm{<mo>/</mo><mo\; stretchy="false">(</mo>}RT\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}{\nu }_0\exp \mathrm{<mo\; stretchy="false">[</mo>}b\mathrm{<mo\; stretchy="false">(</mo>}{T}_0\mathrm{<mo>/</mo>}T-\mathrm{1<mo\; stretchy="false">)</mo><mo>/</mo><mo\; stretchy="false">(</mo>}{T}_0\mathrm{<mo>/</mo>}{T}_u-\mathrm{1<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>},$               (24)
where$T_0\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{T}_u+T_{\ell }\mathrm{<mo\; stretchy="false">)</mo><mo>/</mo>2<mo>=</mo>0.0005}$ . Hence, we have that the activation energy in the model is $E\mathrm{<mo>=</mo>}b{T}_0$  since the gas constant$R\mathrm{<mo>=</mo>1}$ .  In these plots, we can see that the upper TBL is thicker than the lower thermal boundary layer. As the coefficient $b$  and hence activation energy is increased, the thickness of the upper TBL increases. One can see a series of narrow upwelling plumes and downwelling chutes for all of the values of the activation energy. The heat below the upper lid increases slightly (Figures 4 and 5) as the activation energy increases (cf. For the Reynolds model, the increase in heat below the upper lid was much greater with increasing b (Figures 2 and 3)). In models of the Earth, a weak layer called the athenosphere occurs below the 100km thick cold outer lid termed the lithosphere. We see some evidence for a narrow hotter region below the upper lid in the higher activation energy run, and suggest that this temperature anomaly may cause partial melting leading to a weak athenosphere. We also saw a similar but stronger effect in the Reynolds’ model snapshots at high $b$ -values (Figure 2).

A                                                          B

C                                                          D

Figure 4 Snapshots of temperature in simulations using Arrhenius temperature dependence of viscosity for cases equivalent to the Reynolds model with $b\mathrm{<mo>=</mo>7,}b\mathrm{<mo>=</mo>10,}b\mathrm{<mo>=</mo>12}$  and $b\mathrm{<mo>=</mo>15}$  (A−D).⁴

A                                                          B

C                                                          D

Figure 5  Snapshots of mean temperature profiles with depth for simulations using Arrhenius temperature dependence of viscosity for cases of $b\mathrm{<mo>=</mo>7,}b\mathrm{<mo>=</mo>10,}b\mathrm{<mo>=</mo>12}$ and $b\mathrm{<mo>=</mo>15}$  (top-left – bottom-right).

The results here show no pulsation and are similar to the results obtained with Reynolds temperature dependence of viscosity when this dependence is weak, and are similar also, to results in the literature.²⁹ The reason one does not get pulsation in the Arrhenius case is due to the temperature dependence of viscosity being weaker for a given viscosity increase of the upper lid. Namely, as the Arrhenius dependence rises more sharply than the Reynolds dependence, the rate of change of the viscosity at the mean temperature in the model is weaker in the Arrhenius model than in the Reynolds model for a given viscosity drop of the upper lid.

We presented simulations using the Lattice Boltzmann Method, of thermally driven plume dynamics in a 2D model with a hot base and cold lid where the material in the model has a viscosity that is strongly temperature dependent. We find that in extreme cases of strongly temperature dependent viscosity, that plumes seem to pulsate with bubbles of hot low viscosity material pushing upward through the model. However, when the temperature dependence is not so extreme, one observes more classical plume-like structures pushing upward through the model  that compare well with published results.⁴ In the strongly temperature dependent cases, a relatively hot layer forms underneath the cold lid. This may have relevance in explaining the formation of the athenosphere, the weak layer underneath the lithosphere.

DA Yuen would like to thank National Science Foundations geochemistry and CISE program for support and discussions with Matthew Knepley.

Author declares there is no conflict of interest in publishing the article.

1. Schubert G, Turcotte D, Olsen P. Mantle convection in the earth and planets. UK: Cambridge University Press; 2001. p. 913.
2. Yuen DA, Balachandar S, Hansen U. Modelling mantle convection: a significant challenge in geophysical fluid dynamics. In: Fox PA, Kerr RM, editors. Geophysical and Astrophysical Convection. UK: Taylor & Francis Group; 2000. p. 257–94.
3. Hier Majumder CA, Yuen D, Vincent AP. Four dynamical regimes for a starting plume model. Physics of Fluids. 2004;16(5):1516–31.
4. Zhong S, Yuen DA, Moresi LN, et al. Numerical models for mantle convection. Treastise on Geophysics. 2015;7:197–222.
5. Gerya Taras. Introduction to Numerical Geodynamic Modelling. UK: Cambridge University Press; 2009. p. 358.
6. Lallemand P. Hybrid finite-difference thermal lattice Boltzmann equation. Int J Mod Phys B. 2003;17(1):41–7.
7. Schmalzl J, Breuer M, Hansen U. The influence of Prandtl number on the style of vigorous thermal convection. Geophysical & Astrophysical Fluid Dynamics. 2002;40:1–23.
8. Chen S, Doolen G. Lattice Boltzmann method for fluid flows. Annual Review of Fluid Mechanics. 1998;30:329–64.
9. Chapman S, Cowling TG. The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction, and diffusion in gases. UK: Cambridge University Press; 1970. p. 423.
10. Frisch U, Hasslacher B, Pomeau Y. Lattice-gas automata for the Navier-Stokes equation. Phys Rev Lett. 1986;56:1505–8.
11. Bhatnagar PL, Gross EP, Krook M. A model for collision processes in gases. I. Smalled amplitude processes in charged and neutral one-component systems. Phys Rev. 1954;94(3):511–25.
12. Qian Y, D'Humières D, Lallemand P. Lattice BGK models for Navier-Stokes equation. EPL. 1992;17(6):479.
13. Huang H, Sukop MC, Lu XY. Multiphase Lattice Boltzmann Methods: theory and application. USA: Wiley Blackwell; 2015. p. 388.
14. He X, Chen S, Doolen GD. A novel thermal model for the Lattice Boltzmann Method in incompressible limit. J Comp Phys. 1998;146(1):282–300.
15. Guo Z, Shi B, Zheng C. A coupled lattice BGK model for the Boussinesq equations. International Journal for Numerical Methods in Fluids. 2002;39(4):325–42.
16. Shan X. Simulation of Rayleigh-Bénard convection using a Lattice Boltzmann Method. Phys Rev E. 1997;55:2780–8.
17. Wang J, Wang D, Lallemand P, et al. Lattice Boltzmann simulations of thermal convective flows in two dimensions. Computers and Mathematics with Applications. 2013;65(2):262–86.
18. Arun S, Satheesh A, Mohan CG, et al. A review on natural convection heat transfer problems by Lattice Boltzmann Method. J Chem and Pharm Sci. 2017;10(1):635–45.
19. Huber C, Parmigiani A, Chopard B, et al. Lattice Boltzmann model for melting with natural convection. Int J of Heat and Fluid Flow. 2008;29(5):1469–80.
20. Barrios G, Rechtman R, Rojas J. The lattice Boltzmann equation for natural convection in a two-dimensional cavity with a partially heated wall. J Fluid Mech. 2005;522:91–100.
21. Mora P, Yuen DA. Simulation of plume dynamics by the Lattice Boltzmann Method, Geophys. J Int Express Letters. 2017;210(3):1932–7.
22. Bunge HP, Richards MA, Baumgardner JR. Effect of depth-dependent viscosity on the planform of mantle convection. Nature. 1996;379:436–8.
23. Bunge HP, Richards MA, Baumgardner JR. A sensitivity study of 3-dimensional spherical mantle convection at 108 Rayleigh number: Effects of depth-dependent viscosity, heating mode, and an endothermic phase change. J Geophys Res. 1997;102(B6):11991–12007.
24. Wang S, Zhang S, Yuen DA. Visualization of downwellings in 3-D spherical mantle convection. Phys of the Earth and Planetary Interiors. 2007;163(1–4);299–304.
25. Zhang S, Yuen DA. The influences of lower-mantle viscosity stratification on 3-D spherical-shell mantle convection. Earth Planet Sci Lett. 1995;132(1–4):157–66.
26. Shan X, Chen H. Lattice Boltzmann model for simulating flows with multiple phases and components. Phys Rev E. 1993;47(3):1815–9.
27. Mora P, Yuen DA. Simulation of regimes of convection and plume dynamics by the Lattice Boltzmann Method. Physics of the Earth and Planetary Interiors. 2018;275:69–79.
28. Frank Kamenetskii D A. Diffusion and heat transfer in chemical kinetics. USA: Princeton University Press; 1955. p. 382.
29. Zhong S, Zuber MT, Moresi L, et al. Role of temperature dependent viscosity and surface plates in spherical shell models of mantle convection. J Geophys Res. 2000;105(B5):11063–82.