Numerical simulation of mantle convection using a 1 temperature dependent nonlinear viscoelastic model 2 3

8 In the present article, the mantle convection is simulated numerically using a 9 temperature dependent non-linear viscoelastic model for the first time. The numerical 10 domain of problem is considered as a 4000km*2000km rectangular box and the CFD 11 simulation is performed using finite volume method. Unlike the previous works which 12 had been investigated the mantle convection using the linear viscoelastic models or 13 simple nonlinear inelastic viscous equations (such as power law or cross equations), it is 14 solved via the nonlinear Giesekus constitutive equation. Because of large-scale creeping 15 flow in geometry and time, it is shown that the results of Giesekus equation are more 16 reliable for this problem. The main innovative aspects of current study is investigation 17 of temperature dependency of rheological properties of mantle including viscosity, 18 normal stress differences and relaxation time using appropriate equations of state. The 19 variation of gravitational acceleration with depth of Earth and the effect of the work of 20 stress field (viscous dissipation) on mantle convection are also simulated for the first 21 time. 22

that carries out various tectonic activities and continental drift (Bénard (1900), Batchelor (1954), Elder (1968)).This motion occurs on a large scale of space and time.
From fluid mechanics point of view, mantle convection is approximately a known phenomenon; the only force which causes convective flow is buoyancy force while this phenomenon is affected by the nature of non-Newtonian rheology (Christensen (1985)) and depth-and temperature-dependent viscosity.Gurnis and Davies (1986) just used a depth dependent viscosity and assumed that the Rayleigh number is constant.They deduced this phenomenon depend on Rayleigh number, as when Ra is increased, the thermal boundary layer will be thinned and the center of circulation shifts more to the narrow descending limb.Hansen et al. (1993) examined the influences of both depthdependent viscosity and depth-dependent thermal expansivity on the structure of mantle convection using two-dimensional finite-element simulations.They concluded depthdependent properties encourage the formation of a stronger mean flow in the upper mantle, which may be important for promoting long-term polar motions.The rheology of mantle strongly depends on the temperature and hydrostatic pressure (Ranalli (1995), Karato (1997)).Also, because of huge geometry of Earth's mantle (2000km), the gravity cannot be considered as a constant, and it is a function of depth.Kellogg and King (1997) developed a finite element model of convection in a spherical geometry with a temperature-dependent viscosity.They have focused on three different viscosity laws: (1) constant viscosity, (2) weakly temperature-dependent viscosity and (3) strongly temperature-dependent viscosity.Moresi and Solomatov (1995) have simulated it as two-dimensional square cell with free-slip boundaries.They reached an asymptotic regime in the limit of large viscosity contrasts and obtained scaling relations that found to be agreement with theoretical predictions.Ghias and Jarvis (2008) investigated the effects of temperature-and depth-dependent thermal expansivity in two-dimensional mantle convection models.They found the depth and temperature dependence of thermal expansivity each have a significant, but opposite, effect on the mean surface heat flux and the mean surface velocity of the convective system.The effect of temperature-dependent viscosity was studied in literature in twodimensional rectangular domains (Severin and Herwig (1999), Pla et al. (2009), Hirayama and Takaki (1993), Fröhlich et al. (1992)).Tomohiko et al.

GOVERNING EQUATIONS
The governing equations of an incompressible viscoelastic fluid flow consist of the continuity, momentum and energy equations: where U is the velocity vector,  is density, c is heat capacity, p is static pressure, T is temperature, k is thermal conductivity, t is time, u is power of heat source and τ is the total stress tensor.The stress tensor is consisted as the summation of Newtonian n τ and viscoelastic contributions v τ as follows: Since the present study examines mantle convection, this parameter must be near unity.
In other words, the viscoelastic portion dominates to pure Newtonian portion in behavior of fluid flow.Therefore, the main portion of viscosity of mantle could be attributed to the v  .
The Giesekus model is a popular choice, because of its relative success in several flows, and its reduction to several well-known simpler models, which make it useful in a variety of flow situations.The key characteristic of this model is that it includes nonlinear term in stress.Here, the Giesekus model is used as the non-linear constitutive equation: where v  is the viscosity contribution of viscoelastic material at zero shear rate and   1 v τ is the upper convected derivative of viscoelastic stress tensor defined by: . In this paper, the viscosity is assumed to be depended on depth and temperature as follow: where 0  is the total viscosity at reference temperature ( 0 T ), y is the depth (per 1000Km), and  is the exponential rate.The relaxation time (  ) is also assumed to be an exponential function of temperature: Because of large scale of geometry and the nature of mantle convection, the dependency of density on temperature and pressure are considered as follows:

NON-DIMENSIONALIZATION
According to Fig. 1, the Cartesian coordinate system is used in this study.The dimensionless parameters of flow field are as follows: where x and y are indicating the coordinate directions; H is the depth of geometry, W 0 is the reference velocity, 0 is the fluid viscosity,  is density and Re, We and En are the Reynolds, Weissenberg and Elastic numbers, respectively.The ~ notation signifies that parameter has dimension.The governing dimensionless parameters of heat transfer are as follows: .
where  is the thermal expansion coefficient.In order to get closer to reality, in the energy equation, we assume a viscosity dissipation term ( :  τ U ).This term is the  (2002)).In Eq. ( 11c), Φ is the dimensionless form of work of stress field and obtain from following equation: This variation in viscosity introduces a relativity factor in the problem.Here, the nondimensionalization is performed regarding to the value of the viscosity in the upper plate.Therefore, a new Rayleigh number should be defined, due to the variation of viscosity: In our numerical calculations, the values of the parameters are related to the values in the mantle (Pla et al., 2010), Table 1 shows the values of parameters used in calculations.Due to the nature of mantle convection the Pr number and viscosity are assumed to be in order of 10 26 and 20 10 , respectively.Also, a Rayleigh number equal to 227 is used for this simulation.
Remember that the gravitational acceleration of the Earth is decreased by increasing the depth.Because of the large scale of geometry, the variation of gravitational acceleration with depth is considered in present study.For this purpose, we used the data of Bullen (1939) and fitted the following six order interpolation on them with 95% confidence:   where (1000 ) y Km is the depth from bottom plate.We used the above equation in CFD simulation of mantle convection which is the other innovative aspect of present study.

NUMERICAL METHOD, BOUNDARY AND INITIAL CONDITIONS
There are totally eight solution variable parameters in the discretized domains, comprising two velocities and three stress components, pressure, pressure correction and temperature.All of flow parameters are discretized using central differences, except for the convective terms which are approximated by the linear-upwind differencing scheme (LUDS) (Patankar and Spalding (1972)).This is the generalization of the wellknown up-wind differencing scheme (UDS), where the value of a convected variable at a cell face location is given by its value at the first upstream cell center.In the linearupwind differencing scheme, the value of that convected variable at the same cell face is given by a linear extrapolation based on the values of the variable at the two upstream cells.It is, in general, the second-order accurate, as compared with first-order accuracy of UDS, and thus, its use reduces the problem of numerical diffusion (Oliveira et al.
( bottom and top plates are assumed a constant temperature so that the bottom plate has a higher temperature.These boundaries have a zero gradient velocity and tensor components, too.

Grid Study and Validation
We perform some CFD simulations with different number of grids to study the dependency of solution to mesh size.The meshes included quadratic elements. .The numerical error decreases with increasing the number of meshes as the mean error beings less than 0.08% for mesh size greater than140 70  .This finding indicates that a grid-independent solution is obtained when using a mesh sizes larger than 140 70  .To ensure that the obtained solution is grid- independent, a mesh size of 150 75  was used for the CFD simulations.
As a benchmark comparison, simulations for free convection of Newtonian fluid flow between two parallel plate have been carried out at

CFD Simulation of Mantle Convection Using Giesekus Model
In this section, the effects of various parameters on flow regime of mantle convection are studied.As observed in Eq. ( 4), the variation of parameters  and  could affect the stress tensor field and this change in stresses will affect the velocity field.
According to the study of Pla et al. (2010), it could be inferred that with increasing the exponential rate  , the circulations created by natural convection are moved toward the bottom plate.It is resulted from the fact that by increasing  , the viscosity near bottom plate would be decreased and the flow tends to circulate in this place.Also, another parameter that effect on the flow and the circulation intensity is G  .The results of variations of these parameters will discuss in next sections.Remember that the dependency of rheological and thermal properties and density on temperature and pressure are considered and the variation of gravitational acceleration with depth of Earth is modeled in following results.Christensen (1983), Cserepes (1982), Sherburn (2011), Van der Berg (1995), Yoshida (2012)) at n=3, and the Newtonian model used by Pla et al. (2010).This Figure is presented in order to compare the results of current CFD simulation (based on the non-linear Giesekus consecutive equation, thermal-pressure dependence properties and depth dependence gravitational acceleration) with previous simpler simulations that used Newtonian and power-law models.As it is obvious, the velocity near upper plate for Giesekus model is less than from the results of Pla et al. (2010) and power-law model.That is due to the elastic force and higher value of viscosity at lower shear rates.Also, the maximum vertical velocity of our simulation is smaller and the location of maximum vertical velocity occurred upper than the location reported by Pla et al. (2010).That is because of the viscoelastic portion of fluid behavior that we will discuss it in next sections.As it is shown in Fig. 3, the depth in which the maximum velocity occurs is approximately similar for power-law model and Giesekus constitutive equation.That is because of the effect of apparent viscosity dependency to velocity gradient.Also noting to the velocity profile, it is seen that all of models have the same results in vicinity of lower plate.But

Investigation of the Effect of Exponential Rate of Viscosity ( Γ )
We studied firstly the effect of increasing  from zero to 3 10  on mantle convection.This parameter represents the dependency of viscosity on temperature variation.Also, with increasing the value of Γ by 3 10  , it growths up to 2.56 times.Actually, with increasing the exponential rate, the dependency of viscosity on temperature is intensified and then the right hand side of Eq. ( 4) increases so this change leads to enhancement of stress field. .As expected from Eq. ( 6), the viscosity will be more depended on temperature by increasing the value of Γ .Thus, the viscosity will be decreased with increasing Γ and in the other hand; the velocity field will be intensified that the participation of these factors determines stresses in vicinity of upper plate.According to Fig. 9, in the case of 5 Γ 10   , with increasing G  from 0.5 to 0.8, the maximum stress magnitude is increased by 32.2% and by enhancing G  to 0.9 and 0.98, the growing percentages are 32.2% and 101%, respectively.As mentioned before, there are several factors that affect the flow pattern such as Γ and G  .The result of this participation clearly is seen here, when the viscosity ratio vary from 0.9 to 0.98, it seems that in this interval, the effect of these two parameters ( Γ and G  ) is neutralized each other and lead to having the same stress magnitude at these points.

Investigation of the Effect of Viscosity Ratio ( G β )
The parameter G  is a criterion portion for demonstration of domination of viscoelastic towards pure Newtonian portions of fluid behavior.In fact, when this parameter is much closer to unity, the viscoelastic behavior is dominated and when vertical velocity near to the both lower and upper plates is decreased.This effect is related to the higher value of viscosity of viscoelastic potion in comparison of pure Newtonian behavior that causes increasing the total viscosity and decreasing the fluidity of model (refer to Eq. 3).This finding is approved by the data of maximum magnitude of shear stress near to the upper plate which is reported in Table 3.According to the  ).As it is understood from Fig. 11, in constant viscosity ratio, when Γ is increased, the velocities are increasing very strongly, but as viscosity ratio changes, a contrast occurred between these two factors (as it is shown in Fig. 11c, the velocities are increased and in Fig. 11b, the vertical velocities are decreased).In other word, at 0.9 G   , the effect of exponential rate is prevailed but with increasing the viscosity ratio to 0.98 G   , the effect of viscosity ratio is dominated.

Investigation of the Effect of Elasticity
The elastic number is generally used to study the elastic effect on the flow of viscoelastic fluids.According to the Eq. 9, the elastic number is defined as the ratio of Weissenberg to Reynolds numbers.This dimensionless group is independent from kinematic of flow field and it is only depended on material modules for a given  velocity could be attributed to increasing the normal stresses resulted from fluid elasticity.In the other word, some main portion of energy of convection is stored as the elastic normal stresses.In larger elastic numbers, the effective viscosity of flow is decreased which is related to the nature of nonlinear dependency of viscometric function of Giesekus constitutive equation on relaxation time at large enough elastic numbers (Bird et al. (1987))., the magnitudes of normalized velocities in vicinity of upper plate are increasing by enhancing G  from 0.5 to 0.9 between 20% to 50% and with increasing the viscosity ratio to 0.98, the velocities are decreasing about 70%.In contrast, for the lower plate, this variation is reversing, i.e., the velocities with increasing G  to 0.9 are decreasing.The same effect is available for 0.2  

Investigation of Mobility Factor Effect
. Also, the variation of velocity near upper plate for 0.1   and 0.3 are similar.In these cases, with increasing G  from 0.5 to 0.9, the velocities in this place are decreasing and with increasing the viscosity ratio to 0.98, the magnitudes of velocities are ascending.Table 5 presents the maximum normalized vertical velocity for various values of elastic numbers and different viscosity ratios.

Investigation of the Effect of Rayleigh Number
If we want to study natural convection and investigate the strength of convection, the Rayleigh number is a suitable criterion for this aim.Since mantle convection has a low Rayleigh number, thus the temperature field should have a conductive form (see Fig 7).
According to Eq. ( 10), the Rayleigh number is a function of temperature, so it is varying all over the geometry because the viscosity is temperature dependent and is varying.1135, this growth in velocities is in order of 10 2 and when we set the Ra as 2270, the velocity magnitude will be in order of 3 10 .It is important to remember that the temperature difference between the hot and cold plates is the potential of mantle convection so the velocity is increased by increasing the Rayleigh number.Fig. 16 shows the stress contours for various Rayleigh number.The Figure shows that with increasing the Rayleigh number, the maximum stress in geometry has enhanced significantly.This effect is related to increasing the shear rate of flow field which is intensifying the stress field.According to the Figure, the Giesekus model predicts a large shear stress in comparison of normal stress components which is related to the shear flow behavior of mantle convection which has a suitable agreement with previous reports that used other constitutive equations (Ghias and Jarvis (2008), Severin and Herwig (1999), Pla et al. (2009), Hirayama and Takaki (1993), Fröhlich et al. (1992), Tomohiko et al. (2004)).

CONCLUSIONS
Current study deals with a numerical simulation of mantle convection using a temperature dependent nonlinear viscoelastic constitutive equation.The effect of temperature on rheological properties consisting of the viscosity, normal stress It is important to remember that the non-linear constitutive equations such as the Giesekus equation could able to model the material elasticity and relaxation spectra much better than linear models for large deformations of flow field.We also showed that the result of this model has an obvious deviation from pure Newtonian and powerlaw solutions that reported in literatures.
The effect of temperature on viscosity of the mantle is studied, firstly.The results show that increasing of exponential viscosity rate led to the enhancing the maximum velocity and making the circulation moving downward so that with increasing  from zero to 10 -3 , an increase of 4.32 times in vertical velocity and an increase of 2.56 times in xx  were obtained.A formula have presented for the position of maximum vertical velocity as a function of  .The effect of viscosity ratio is also investigated on the In constant viscosity ratio, when G  increases, the velocities are rising very strongly, but as viscosity ratio changes, a competition occurred between these two factors.In other word, at 0.9 G   , the effect of exponential rate is prevailed but with increasing Table 1.Parameters related to mantle convection (Pla et al. (2010)).
a creeping flow in the mantle of the Earth that causes some 27 convective currents in it and transfers heat between core and Earth's surface.In fluid 28 mechanics, the free convection is a classic topic driven by the effect of temperature 29 gradient on density.This solid-state convection in mantle is an abstruse phenomenon 30 Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.
(2004)  simulated a two-dimensional rectangular domain with assuming the mantle as an incompressible fluid with a power-law viscosity model.They employed a simplified two-layer conductivity model and studied the effects of depth-dependent thermal conductivity on convection using two-dimensional Boussinesq convection model with an infinite Prandtl number.Their results implied that the particular values of thermal conductivity in horizontal boundaries could exert more significant influence on convection than the thermal conductivity in the mantle interior.Stein et al. (2004) explored the effect of different aspect ratios and a stress-and pressure-dependent viscosity on mantle convection using three-dimensional numerical simulation.Ozbench et al. (2008) presented a model of large-scale mantle-lithosphere dynamics with a temperaturedependent viscosity.Ichikawa et al. (2013) simulated a time-dependent convection of fluid under the extended Boussinesq approximation in a model of two-dimensional rectangular box with a temperature-and pressure-dependent viscosity and a viscoplastic property.Stien and Hansen (2008) employed a three-dimensional mantle convection Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.model with a strong temperature, pressure and stress dependence of viscosity and they used a viscoplastic rheology.Kameyama and Ogawa (2000) solved thermal convection of a Newtonian fluid with temperature-dependent viscosity in a two-dimensional rectangular box.Kameyama et al. (2008) considered a thermal convection of a high viscous and incompressible fluid with a variable Newtonian viscosity in a threedimensional spherical geometry.Gerya and Yuen (2007) simulated a two-dimensional geometry and non-Newtonian rheology using power-law model.In the present paper, the mantle convection is simulated numerically using a temperature dependent non-linear viscoelastic model for the first time.The geometry of problem is shown in Fig. 1.Here, the calculation domain is considered as a 4000km×2000km rectangular box.Two hot and cold plates are considered at the bottom and top of box, respectively.The isolator thermal condition is considered at the left and right hand sides of domain.The problem is solved via a second order finite volume method.The effect of temperature on rheological properties consist of the viscosity, normal stress differences and relaxation time of mantle are modeled using appropriate equations of state which are the main innovative aspects of current study.The variation of gravitational acceleration with depth of Earth and the effect of the work of stress field (viscous dissipation) on mantle convection are simulated for the first time.According to the literature, the previous studies are restricted to the linear and quasi-linear viscoelastic constitutive equations and the nonlinearity nature of mantle convection was modeled as simple nonlinear constitutive equations just for apparent viscosity such as the power-law and cross models.Here, the Giesekus nonlinear viscoelastic model is Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.used as the constitutive equation.

)
Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License. is derived by kinetic theory, arising naturally for polymer solutions.This model contains four parameters: a relaxation time  ; the solvent and polymeric contributions at the zero-shear rate viscosity, n  and v  ; and the dimensionless "mobility factor"  (Bird et al. (1987)).The origin of the term involving  can be associated with anisotropic Brownian motion and/or anisotropic hydrodynamic drag on the constitutive of heavy particles.


is density at reference temperature and pressure,  is thermal expansivity and C  is compressibility coefficient.
10) Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.temperature of fluid, respectively; k is the conduction coefficient,  is thermal diffusivity, h is the convection heat transfer coefficient and Br , Pr , Ra and Nu are the Brinkman, Prandtl, Rayleigh and Nusselt numbers, respectively.Thus, the dimensionless form of continuity and momentum equations are as follows: effect of stress field work on fluid flow and for Newtonian fluids; it has always a positive sign according to the second law of thermodynamic.Actually, this positive term refer to the irreversibility of flow field work and thus in Newtonian fluid it is known as viscosity dissipation.The interesting point of this term for viscoelastic fluids is the local possibility of being negative.In effect, having locally negative value of this term indicates that part of energy is saved in elastic constituent of fluid (Bird et al.Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.
)).The Cartesian reference coordinate system is located in the bottom boundary and at left corner.Boundary conditions consist of two adiabatic walls in west and east and two isothermal walls at north and south.For all boundaries, a no-slip condition is imposed for the fluid velocity.The rest situation is used as the initial condition.The used geometry and boundary conditions in this study are shown in Fig.1.The geometry has a rectangular shape with an aspect ratio of 2. Boundary conditions consist of two isolated walls with zero gradient stress tensor components.The boundary conditions for Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.
previously by Khezar et al. (2012) and Turan et al. (2011) for power-law fluid.The diagrams of average Nusselt number obtained from the present study and work of Khezar et al. (2012) at n=1 are shown in Fig. 2a.As an additional benchmark comparison, the distribution of dimensionless vertical velocity reported by Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.Turan et al. (2011) and the results obtained from the present study are illustrated in Fig.It is understood that in both cases, the results of present CFD simulation have a suitable agreement with results of Khezar et al.(2012) andTuran et al. (2011) with maximum error less than 3%.

Solid
Fig.3demonstrates a comparison between vertical velocity profiles of our nonlinear viscoelastic model, power-law model (reported byChristensen (1983), for upper plate, the Figure demonstrates that the slope of vertical velocity for the Giesekus model is smaller than the others.According to the Figure, there is a resistance against the upward flow for Giesekus profile that two other models cannot predict it.Actually, that is due to the consideration of elastic portion of fluid flow in our numerical simulation.This finding indicated that the velocity and stress field have an obvious deviation from Newtonian and generalized Newtonian behaviors by considering a non-Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.linear constitutive equation for mantle convection.In next sections, the effects of material and thermal modules on mantle convection are studied based on the CFD simulations that obtained using Giesekus non-linear model.


evident from Fig.4that the circulations in the mantle physically depend on  .As the exponential rate ( Γ ) is increased, the maximum velocity in geometry is enhanced and the circulations moved downward.According to Eq. (6), the dependency of viscosity of mantle on temperature is more increased by enhancing the exponential rate ( Γ ).In other words, by increasing the exponential rate ( Γ ), the viscosity is more decreased near to the lower plate (high temperature region) and the fluency of mantle is intensified.Therefore, it is expected that the velocity of mantle convection is enhanced by increasing the exponential rate.The results show that an increment of 1.6% in vertical velocities by increasing the exponential rate from zero to 5 10  , 17.1% growth by increasing Γ to 4 10  and with enhancing the  from zero to 3 10  it growths up to 4.32 times.The CFD simulations indicated that the effect of exponential rate on maximum value of velocity is nonlinear.The contours of axial normal stress and shear stress are shown in Fig.5.As it is obvious, the exponential rate Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.has a significant influence on magnitude of stress fields that is increased by enhancing the exponential rate.As an example, becomes 1.1 times greater than the one with exponential rate of zero.

Fig
The above correlation is used in plotting the Fig.6.The downward movement of location of maximum vertical velocity with increasing exponential rate could be attributed to shifting the center of vortices which is shown previously in Fig.5.In Fig.7, the temperature distribution in mantle is shown.According to this Figure, heat transfer regime is almost conduction.Nevertheless, closer looking to the temperature distribution, some convection behavior could be observed.The temperature profile on a horizontal line is shown in Fig. 8.As it is expected, the temperature profile shown in Fig. 8 has a minimum value at mid of horizontal line and the maximum values Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.are located at left and right hand sides of numerical domain.Fig. 9 shows the stress magnitude on upper plate for different value of Γ

G
is close to zero, the pure Newtonian behavior of fluid is dominated.As it is shown in Fig. 10, by increasing G  from 0.8 to 0.98, the stress magnitude on upper plate has been increased, but the Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.

Fig. 11
Fig. 11 shows variation of normalized vertical velocity on a vertical line for different values of exponential rates ( Γ ) and viscosity ratios (

Solid
Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.geometry.Here, the elastic number is proportional with relaxation time of model and it is increased by enhancing the material elasticity.Figs. 12 and 13 display velocity and stress magnitude for different values of elastic number.

Fig. 14
Fig.14shows the effects of mobility factor on the vertical velocity for different values of viscosity ratio.Due to the non-linear nature of our viscoelastic model and the high elastic number, anticipation of effects of all factors is not easy and it is strongly affected by the variation of other factors.Regarding to high viscosity of mantle, the effect of mobility factor must be minimal, as it is shown in Fig.14.The effects of mobility factor are only important near both upper and lower plate.In the other word, the main variation of velocity distributions with changing the mobility factor occurs in the upper

Fig. 15
Fig.15 presents the streamlines for different Rayleigh numbers.According to Fig. 15, by increasing the Rayleigh number, the velocity in geometry is increased and the circulations move downward and get more intense.By increasing Ra from 22.7 to 227, the velocity magnitude will vary with order of 10 1 .If we rise the Rayleigh number to differences and relaxation time of mantle are modeled using appropriate equations of state which were the main innovative aspects of current study.The variation of gravitational acceleration with depth of Earth and the effect of the work of stress field Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.(viscous dissipation) on mantle convection were simulated for the first time.According to the literature, the previous studies were restricted to the linear and quasi-linear viscoelastic constitutive equations and the nonlinearity nature of mantle convection was modeled using simple nonlinear constitutive equations just for apparent viscosity such as the power-law and cross models.The Giesekus nonlinear viscoelastic model was used as the constitutive equation in present study.This high order nonlinear model was used because of large-scale creeping viscoelastic flow of mantle convection in space and time.Using Giesekus constitutive equation, we present a more accurate solution for this problem because of taking into account of shear-dependent nonlinear viscometric functions, the effects of third invariant of shear rate tensor on stress field, and effects of material elasticity for large deformations of mantle.

Solid
Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.mantle convection.These results not only show how stress magnitude on upper plate increases by enhancing the viscosity ratio from 0.8 to 0.98, but also prove decreasing of the vertical velocity near to the both lower and upper plates.These effects are related to the higher value of viscosity of viscoelastic Gesikus model relative to the pure viscous portion (Newtonian behavior) which causes decreasing of fluidity of mantle convection.
viscosity ratio is dominated and the velocities are descended.The variation of Elastic number shows the nature of nonlinear dependency of viscometric function of Giesekus constitutive equations on relaxation time at large enough elastic numbers.Present study indicates decreasing of effective viscosity flow for larger elastic numbers.The obtained results show how main variations of velocity distributions with changing of mobility factor occur in the upper and lower plates.Here, the effect of Rayleigh number on mantle convection is also investigated and characterized that with increasing the Rayleigh number, the maximum stress in geometry has enhanced significantly.This effect is related to increasing the shear rate of flow field which is intensifying the stress field.Future works could be focused on the effect of mantle convection on plate motions, effect of chemical reactions occurring in the mantle, and plumes growing by considering a non-linear viscoelastic consecutive equation.Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.573 574 ., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.
This high order nonlinear model is used because of large-scale creeping viscoelastic flow of mantle convection in domain and time.Using

Table , max
is increased by enhancing the viscosity ratio which is caused from increasing the fluid viscosity.
Table 4 presents the value of maximum normalized vertical velocity for different elastic numbers and various viscosity ratios.According to the Fig. 12, the velocity of mantle convection is decreased

Table 2 .
Percentage of mean absolute errors between average velocity obtained from different meshes and the 200 100  reference mesh.Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.

Table 3 .
Maximum magnitude of stress on top plate for different values of G Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License. G

Table 4 .
Maximum magnitude of vertical velocity on a vertical line at x=1 for different Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License. G

Table 5 .
Maximum magnitude of vertical velocity on a vertical line at x=1 for different Solid Earth Discuss., doi:10.5194/se-2016-12,2016 Manuscript under review for journal Solid Earth Published: 15 February 2016 c Author(s) 2016.CC-BY 3.0 License.