Mathematical modelling of mantle convection at a high Rayleigh number with variable viscosity and viscous dissipation

Introduction Convection in mantle is responsible for most of the physical and chemical phenomena happening on the surface and in the interior of the Earth, and it is caused by the heat transfer from the interior to the Earth’s surface. Even though there are some debates, it is quite well established that convection in the mantle is the driving mechanism for plate tectonics, seafloor spreading, volcanic eruptions, earthquakes, etc. [1]. However, the mechanism of mantle convection is still an unsolved mystery since the rheology of mantle rocks is extremely complicated [2–4]. Temperature, pressure, stress, radiogenic elements, creep, and many other factors influence the mantle’s behavior on a large scale. One of its significant but complex characteristics is its viscosity, which is dependent mainly on temperature, pressure, and stress [5]. In earlier studies of mantle convection, scientists assumed constant viscosity (e.g. [6, 7]) but later, among many others Moresi and Solomatov [8, 9], studied the temperature-dependent viscosity case numerically and concluded that the formation of an immobile lithosphere on terrestrial planets like Mars and Venus seems to be a natural result of temperature-dependent viscosity. However, studies with purely temperature-dependent viscosity cannot portray the true convection Abstract

pattern of the Earth's mantle. As a result, convection with temperature and pressuredependent viscosity is becoming more important, and some notable works in this area have recently been published [10][11][12][13][14]. Christensen [10] showed that additional pressure dependence of viscosity strongly influences the flow regimes. In a 2D axi-symmetrical model, Shahraki and Schmeling [15] examined the simultaneous effect of pressure and temperature-dependent rheology on convection and geoid above the plumes, and Fowler et al. [16] studied the asymptotic structure of mantle convection at high viscosity contrast.
According to King et al. [17], when pressure increases through the mantle, there is a corresponding increase in density due to self-compression. In a vigorously convecting mantle, the rate at which viscous dissipation, which is the irreversible process that changes other forces into heat, is non-negligible and contributes to the heat energy of the fluid, resulting in adiabatic temperature and density gradients that reduce the vigour of convection. Conrad and Hager [18] proposed that the viscous dissipation and resisting force to plate motion may have significant effects on convection and the thermal evolution history of the Earth's mantle. Leng and Zhong [19] concluded that the dissipation occurring in a subduction zone is 10-20% of the total dissipation for cases with only temperature-dependent viscosity, whereas Morgan et al. [20] declared that when slabs subduct, about 86% of the gravitational energy for the whole mantle flow is mostly transformed into heat by viscous dissipation. According to Balachandar et al. [21], numerical simulations of 3D convection with temperature-dependent viscosity and viscous heating at realistic Rayleigh numbers for Earth's mantle reveal that, in the strongly timedependent regime, very intense localized heating takes place along the top portion of descending cold sheets and also at locations where the ascending plume heads impinge at the surface. They also found that the horizontally averaged viscous dissipation is concentrated at the top of the convecting layer and has a magnitude comparable to that of radioactive heating. King et al. [17] worked on a benchmark for 2-D Cartesian compressible convection in the Earth's mantle where they used steady-state constant and temperature-dependent viscosity cases as well as time-dependent constant viscosity cases. In their work, the Rayleigh numbers are near 10 6 and dissipation numbers are between 0 and 2, and they conclude that the most unstable wavelengths of compressible convection are smaller than those of incompressible convection. As the research on mantle convection is growing, the importance of studying viscous dissipation is also increasing since it was suggested that the bending of long and highly viscous plates at subduction zones dissipates most of the energy that drives mantle convection [22]. Some notable recent works on numerical studies of convection and effects of variable viscosity and viscous dissipation have been done by Ushachew et al. [23], Megahed [24], Ferdows et al. [25], Ahmed et al. [26], Fetecau et al. [27].
Although mantle convection is a 3D problem, many 2D codes have been developed to gain an understanding of the fundamental mechanism and to minimize the computational cost and complexity. As the Earth's mantle has been affected by many complexities, its basic understanding has been constructed through research on simple Rayleigh-Bénard convection [2]. Over the years, the Rayleigh-Bénard convection has become a benchmark problem in computational geophysics as a paradigm for convection in the Earth's mantle. Although Rayleigh-Benard convection with viscosity variation is a well-known topic for mantle convection, very high viscosity variation (up to 10 30 ) for mantle convection is not widely covered. To the best of our knowledge, mantle convection with strongly variable viscosity, which is temperature dependent and also both temperature and pressure dependent with the inclusion of viscous dissipation, has not been studied so far. The governing equation in two-dimensional form ensures the conservation of mass, momentum, and energy and the thermodynamic equation of state. In this study, incompressible mantle convection will be considered where the mantle viscosity depends strongly on both temperature and pressure, and viscous dissipation is also considered. The convection will be investigated at a high Rayleigh number with high viscosity variations across the mantle.
In "Methods" section the full governing equations for mantle convection and the appropriate boundary conditions for classical Rayleigh-Bénard convection in a 2D square cell are described. The equations are non-dimensionalized and the dimensionless parameters are identified. Though the variable viscosity is defined in an Arrhenius form, a modified form of viscosity is used to improve the efficiency of numerical computation. The computational method for simulation is also described, and the code is verified using some benchmark values. Then the governing model is solved numerically in a unit aspect-ratio cell for extremely large viscosity variations, and steady solutions for temperature and streamlines are obtained. The numerical and graphical results of the computation are described in "Result and discussion" section. Finally, in "Conclusion" section some concluding remarks on the results are given.

Governing equations
A classical Rayleigh-Bénard convection in a two-dimensional unit aspect ratio cell with a free slip boundary condition is taken into account. The temperature difference is fixed between the horizontal boundaries. The convective cell is assumed to be a section of a periodic structure in the associated infinite horizontal layer. When adopting Cartesian coordinates (x, z) with horizontal x-axis and vertical z-axis, the Boussinesq approximation is assumed, which suggests that density variation is barely vital within the buoyancy term of the momentum equation, so that mass conservation takes the shape of the incompressibility condition [16]. The inertia terms within the Navier-Stokes equations (taking the limit of an infinite Prandtl number) are neglected as well. According to Solomatov [28], the integral viscous dissipation within the layer is often balanced by the integral mechanical work done by thermal convection per unit time, and if the viscosity contrast is large, dissipation in the cold boundary layer becomes comparable with the dissipation in the internal region. Thus, in order to balance the energy equation, the extended Boussinesq approximation is used. Here, "extended Boussinesq approximation" means that apart from the driving buoyancy forces, the fluid is treated as being incompressible all over. The non-Boussinesq effects of the adiabatic gradient and frictional heating are introduced into the energy equation [29]. The governing equations ensure the conservation of mass, momentum, and energy. This also ensures a suitable thermodynamic equation of state. The Navier-Stokes equations, which describe the motion in component forms, are [30]  The energy equation is Here, P is the pressure, τ is the deviatoric stress tensor, t is time, ρ is the density, u = (u, 0, w) is the fluid velocity, where u and w are velocity components in the x-and z-directions, g is the assumed constant gravitational acceleration acting downwards (the variation of g across the mantle is quite small that it is taken as constant), τ 1 and τ 3 are the longitudinal and shear components of the deviatoric stress tensor, respectively, η is the viscosity, T b is the basal temperature, ρ 0 is the basal density, κ is the thermal diffusivity, T is the absolute temperature, C p is the specific heat at constant pressure, and α is the thermal expansion coefficient. The deviatoric stress tensor, τ can be expressed as where τ 1 and τ 3 are the longitudinal and shear components of the deviatoric stress tensor, respectively.

The Arrhenius form of viscosity function is
where A is the rate factor, n is the flow index, E is the activation energy, V is the activation volume, and R is the universal gas constant [5]. A unit aspect-ratio cell with a free-slip boundary condition is considered. The temperatures at the bottom and top boundaries are taken as constant, and thermal insulation is assumed on the side walls. The boundary conditions are where d is the depth of the convection cell, T b and T s are the basal and top temperatures, respectively ( Fig. 1). Throughout this work, Newtonian rheology is considered with n = 1 in the viscosity relation and internal heating is neglected. To see the effects of variable viscosity (both temperature-dependent and temperature-and pressure-dependent viscosity) and viscous dissipation on convection, these assumptions are made to make the model less complicated.

Non-dimensionalization and simplification
In order to non-dimensionalize the model, the variables are set as [7,30] Using these in equations from (1) to (3) and dropping the asterisk decorations, the dimensionless equations becomes while the dimensionless version of constitutive relation (3) reads  30:5 in which the dimensionless parameters are, Since this model was developed for the mantle, the typical values of the parameters are given in Table 1, and it is found that for Ra >> 1 , B /Ra can be easily ignored. Therefore, the dimensionless energy equation (7) becomes and viscosity relation (8) becomes Gravitational acceleration g 10 m s −2 Temperature at the top of the mantle T s 300 K Temperature at the base of the mantle  (4) become The dimensionless model consists of governing Eqs. (6), (10), viscosity relation (11) and boundary conditions (12).

Low temperature cut-off viscosity
To investigate the convection with extremely high viscosity contrasts in the mantle layer, a low temperature cut-off viscosity function is used. This cut-off viscosity relation helps reduce the computational stiffness while retaining the sensitivity of the viscosity to the changes in temperature and pressure across the mantle. It is a well-established fact that in strongly temperature-dependent viscous convection, most of the viscosity variation occurs in a stagnant lid in which the velocity is essentially zero. Based on this fact, the sub-lid convection field is calculated accurately (but not the stress field) by cutting off the dimensionless viscosity at a sufficiently high value that the lid thickness, which essentially only depends on the interaction of the lid temperature with the underlying convection flow, is unaffected. The low temperature cut-off viscosity function has the following form where and the cut-off viscosity value 10 r is to be chosen appropriately; in numerical experiments, it is chosen r = 6 . Similar type of Arrhenius law with an imposed cut-off viscosity was applied by Huang et al. [31], Huang and Zhong [32], King [33] and Khaleque et al. [13]. A comparison between full-form viscosity function and cut-off viscosity function is shown in "Comparison with benchmark values and validation" section. Three useful diagnostic quantities which will be used to characterize are viscosity contrast, Nusselt number and root mean square velocity respectively.
The viscosity contrast �η is the ratio between the surface and basal values of the viscosity, defined as (11) The Nusselt number Nu is the ratio of the average surface heat flow from the convective solution to the heat flow due to conduction. It is calculated in the present case of a square cell by the dimensionless relation Nu is equal to unity for conduction and exceeds unity as soon as convection starts.
The vigour of the circulating flow is characterised by the non-dimensional RMS (root mean square) velocity. Here RMS velocity is defined by where u is the horizontal component of velocity and w is the vertical component of velocity.

Computational method
In order to solve the dimensionless governing Eqs. where E i is the estimated error and ε = 10 −6 . Further details of the method can be found in Zimmerman [34].

Comparison with benchmark values and validation
The values of Nusselt number Nu and root mean square velocity V rms are compared with the benchmark values from Blankenbach et al. [35] a and Koglin Jr et al. [36] Table 2 for constant viscosity case. Their values were computed for Ra up to 10 6 and 10 7 respectively. From Table 2, it is evident that the agreement is within a very good range. Then the computation is done with variable viscosity with a high viscosity contrast across the mantle layer. The values of Nusselt number Nu that are compared in Table 3 are found using the full form viscosity function (11) and the cut-off viscosity function (13) for µ = 0.5 and µ = 0.0 . It should be noted that µ = 0.0 indicates temperaturedependent viscosity, whereas µ = 0 implies that viscosity depends on both temperature and pressure. From Table 3 it can be seen that the values of Nusselt number, Nu with full form viscosity function and the values of Nusselt number, Nu with cut-off viscosity function are very close, which validates the use of the cut-off viscosity function for numerical computation.
Tables 4 and 5 show that for each fixed value of µ and D, Nu and V rms decrease as the viscosity contrast increases (i.e., the temperature dependence parameter decreases) across the mantle. It confirms that at the higher viscosity variation, convection becomes weaker, which can also be seen clearly in the thermal distribution Figs. 2 and 3. Nu and V rms values also decrease as D increases for every particular value of µ.  Table 3 Comparison of Nusselt number, Nu of full-form viscosity function (11) and cut-off viscosity function (13)   It is also observed that at a specific viscosity contrast as the pressure dependence parameter µ is increased, both Nu and V rms values increase for a fixed dissipation number D = 0.3 and D = 0.6 . The reason behind this is that even though µ is increased, ε is actually decreased to maintain the fixed viscosity contrast. However, for D = 0.0 , the trend is not that smooth at higher viscosity variations.  In panel 2a, b and 3a, b, the viscosity depends only on temperature (i.e. µ=0.0) and in panel 2c, f and 3c, f, the viscosity depends on both temperature and pressure (i.e. µ = 0.0 ). At each plot of the temperature profile, the blue region corresponds to the cooler temperature whereas the red region corresponds to the high temperature.
For µ = 0.0 , µ = 0.5 and µ = 1.0 , from Figs. 2 and 3 we see that as viscosity contrast �η increases the thickness of the cold thermal boundary layer at the top of the cell. At the lower mantle, which is near the core of the Earth, the boundary is hot as the temperature is very high and this temperature continues to increase as the viscosity contrast gets larger.      The interior temperature decreases significantly as the pressure dependence parameter is included. The convection cell is quite different when viscosity is both temperature and pressure dependent rather than only temperature dependent. Compared to µ = 0.5 the significance of pressure can be seen clearly for µ = 1.0 from both Figs. 2 and 3. The stream function contours where stream function �(x, z) defined as are presented in Fig. 4 for D = 0.3 . As the streamlines represent fluid flow, the absence of a streamline confirms that fluid in that region is immobile. In other words, this immobile region represents the stagnant lid. With increasing viscosity contrast and viscous dissipation, the changes in the convection pattern are very clear. It is observed that the cold thermal boundary layer thickness increases with viscosity contrast. But for a fixed dissipation number, the cold thermal boundary thickness is reduced with the inclusion of the pressure-dependent parameter µ . Clearly, the lid thickness decreases as the pressure dependence parameter is increased at a fixed viscosity variation. However, the lid thickness increases when viscosity variation is increased at a fixed pressure dependence parameter µ and dissipation number D. The Tables 4 and 5 clearly indicate that the heat transfer rate and the root mean square velocity decrease, and Figs. 2, 3 and 4 show that the immobile lid thickness increases as the viscosity contrast at a fixed pressure dependent parameter is increased. The decrease in Nu and V rms values, as well as the increase in the thickness of the cold thermal boundary layer, imply that the convection becomes significantly weaker. A visualization of the isothermal contours in Fig. 5 shows that the hot thermal boundary layer is very thin compared to the cold thermal boundary layer. This figure represents the isothermal contours of a convection cell with temperature dependent viscosity at different viscosity contrast (i.e �η = 10 15 and �η = 10 30 ) when viscous dissipation numbers are D = 0.3 and D = 0.5 . There might not be any significant difference in the convection pattern (i.e., isothermal contours), but the contours are not similar. They are clearly affected by different viscous dissipation numbers at different viscosity contrasts. Isothermal contours (Fig. 6a) and viscosity distribution (Fig. 6b) for µ = 1.0 at �η = 10 30 and viscous dissipation D = 0.3 are shown in Fig. 6. The viscosity variation from top to bottom is shown in Fig. 6b, and the resulting color ranges from the lowest value (blue) to 10 6 (brown). Clearly, the cut-off viscosity function simply ignores the high value of the lid viscosity and considers it as a constant there. Figure 6b shows a low viscosity region in the upper mantle and a relatively high viscosity region in the lower mantle just above the bottom boundary layer. This implies that the interior is not isoviscous.
Horizontally average temperature vs depth profiles for viscous dissipation of D = 0.3 and D = 0.6 are presented in Fig. 7. These figures show how the horizontally averaged temperature varies with depth at different viscous dissipation numbers and at different viscosity variations. It also shows how it changes for temperature-dependent viscosity and temperature-and pressure-dependent viscosity. The rapid change in temperature near the cold upper boundary and the hot lower boundary explains the strong temperature gradients in those regions. The plots also indicate that the core of the mantle, i.e. the interior, is not isothermal for both the temperature dependent viscosity case and the temperature and pressure dependent viscosity case. The interior of the convection cell undergoes a larger jump in temperature when dissipation effect is stronger ( D = 0.6 ). The figures show that the interior temperature increases with the increase of viscosity contrast across the mantle layer for µ = 0.0 and µ = 0.5 at D = 0.3 and D = 0.6 . Similar situation occurs for µ = 1.0 at D = 0.6 but when D = 0.3 , temperature decreases at high viscosity contrast (i.e. at �η = 10 30 ).

Conclusion
The study of a basally heated convection model with a strongly temperature and pressure dependent viscous fluid relative to the Earth's mantle in the presence of viscous dissipation has been the principal aim of this work. The classical Rayleigh-Bénard convection model was solved using a low temperature cut-off viscosity function to avoid the stiffness of computation. It was aimed to pursue viscosity that is dependent only on temperature and simultaneously dependent on both temperature and pressure, and a comparison is presented through figures and tables.
According to Jarvis and Mckenzie [37], the dissipation number is between 0.25 and 0.8, whereas Leng and Zhong [19] estimate D to be 0.5 to 0.7. Ricard [38] found that its value is about 1.0 near the surface, and decreases to about 0.2 near the CMB. From Table 1, D ≈ 0.6 has been found. Thus, the effect of various viscous dissipation numbers for mantle like convection with Ra = 10 7 is checked. The different values of viscous dissipation number show the changes in heat transfer rate Nu and root mean square velocity V rms . It is shown that the fluid is not isothermal and isoviscous in the presence of viscous dissipation in both cases when viscosity is temperature-dependent and temperaturepressure-dependent. The viscosity distribution at high viscosity contrast for µ = 1.0 also showed that the fluid is not isoviscous. Analysis of the results can predict that if the dissipation number is increased, the lid thickness will increase more and the convection rate will decrease notably. But it is also clear that the inclusion of viscous dissipation does not affect the convection pattern in any drastic way. The convection becomes weaker as viscosity contrast becomes larger and the viscous dissipation number is increased. However, the variation in Nu, V rms increase as µ goes from 0 to 0.5, but the trend is different when µ goes from 0.5 to 1.0. Thus, strong pressure dependence in viscosity affects the convection in a different way. For a temperature-dependent viscosity case and a temperature and pressure-dependent viscosity case, the horizontally averaged temperature increases with viscosity contrast in the interior, but the trend is opposite in the top boundary layer, i.e., the stagnant lid. In this study we investigated convection with high viscosity contrast, because for the typical parameter values, it is estimated that the viscosity contrast for the Earth's mantle is 10 50 or more. Without extreme parameter values, it is quite impossible to obtain a proper asymptotic structure of mantle convection for the Earth and other planets. Thus, it is believed that this study will have a significant impact on the study of thermal convection in the Earth's mantle and other planets where viscosity is strongly variable and the variation of the order of magnitude is very large.