Bifurcations and dynamics in convection with temperature-dependent viscosity in the presence of the O(2) symmetry
BBifurcations and dynamics in convection withtemperature-dependent viscosity in the presence of theO(2) symmetry
J. Curbelo , , A. M. Mancho Instituto de Ciencias Matem´aticas (CSIC-UAM-UCM-UC3M),Nicol´as Cabrera, 13-15, 28049 Madrid, Spain Departamento de Matem´aticas, Facultad de Ciencias,Universidad Aut´onoma de Madrid, 28049 Madrid, SpainOctober 30, 2018
Abstract
We focus the study of a convection problem in a 2D set–up in the presence of the O(2)symmetry. The viscosity in the fluid depends on the temperature as it changes its valueabruptly in an interval around a temperature of transition. The influence of the viscositylaw on the morphology of the plumes is examined for several parameter settings, anda variety of shapes ranging from spout to mushroom shaped is found. We explore theimpact of the symmetry on the time evolution of this type of fluid, and find solutionswhich are greatly influenced by its presence: at a large aspect ratio and high Rayleighnumbers, traveling waves, heteroclinic connections and chaotic regimes are found. Thesesolutions, which are due to the symmetry presence, have not been previously describedin the context of temperature dependent viscosities. However, similarities are found withsolutions described in other contexts such as flame propagation problems or convectionproblems with constant viscosity also under the presence of the O(2) symmetry, thusconfirming the determining role of the symmetry in the dynamics.
This paper addresses the numerical study of convection at infinite Prandtl number in flu-ids in which viscosity strongly depends on temperature in the presence of O(2) symmetry.Convection in fluids with temperature-dependent viscosity is of interest because of its impor-tance in engineering and geophysics. Linear and quadratic dependencies of the viscosity ontemperature have been discussed [47, 52, 22, 65], but in order to address the Earth’s uppermantle convection, in which viscosity contrasts are of several orders of magnitude, a strongerdependence with temperature must be considered. This problem has been approached bothin experiments [53, 13, 6, 67] and in theory [45, 7, 50, 44, 60, 61]. In these contexts, thedependence of viscosity with temperature is expressed by means of an Arrhenius law. In [7],the exponential dependence is discussed as an approach to the Arrhenius law by means ofa Taylor series around a reference temperature. This is also called the Frank-Kamenetskiiapproximation (see [27]). Viscosity has also been considered when it depends on other mag-nitudes such as depth [10, 5], a combination of both depth and temperature [5] or pressure[51]. However, it is commonly accepted [20, 51] that in the Earth’s interior, viscosity dependsmost significantly on temperature. The usual approach in numerical models of the mantle[44, 5, 63] is to consider constant thermal conductivity. This approach has also been veri-fied in fluid experiments seeking to model mantle convection [6]. However, studies also existwhich consider variations on thermal conductivity [23, 24, 68].1 a r X i v : . [ phy s i c s . f l u - dyn ] N ov ere, we focus on the study of a fluid in which the viscosity changes abruptly in atemperature interval around a temperature of transition. This defines a phase change overa mushy region, which expresses the melting of minerals or other components. Meltingand solidification processes are important in magma chamber dynamics [8, 9], in volcanicconduits [25, 42], in the formation of chimneys in mushy layers [15], in metal processing inindustry (see, for example, [55]), etc. In phase transitions, other fluid properties in addition toviscosity may change abruptly, such as density or thermal diffusivity. However, in this studywe consider solely the study of effects due to the variability of viscosity, since consideration ofthe effect of simultaneous variations on all the properties prevents a focused understandingof the exact role played by each one of these properties. Viscosity is a measure of fluidresistance to gradual deformation, and in this sense very viscous fluids are more likely tobehave rigidly when compared to less viscous fluids. When examining the proposed transitionwith temperature, we focus on the global fluid motion when some parts of it tend to be morerigid than others. Disregarding the variations on density in this transition moves us awayfrom instabilities caused by abrupt density changes such as the Rayleigh Taylor instability,in which a denser fluid over a lighter one tends to penetrate it by forming a fingering pattern.A recent article by M. Ulvrov´a et al. [64] deals with a similar problem to ours, but takes intoaccount both variations in density and in viscosity. Thermal conductivity effects are relatedto the relative importance of heat advection versus diffusion. In this way, diffusive effects areimportant at large conductivity, while heat advection by fluid particles is dominant at lowconductivity. The contrasts arising from these variations are beyond the scope of our workand thus are disregarded here.This paper addresses the convection of a 2D fluid layer with temperature-dependent vis-cosity and periodic boundary conditions possessing the O(2) symmetry . The motivation arisesfrom the fact that symmetric systems typically exhibit more complicated bifurcations thannon-symmetric systems and introduce conditions and degeneracies in bifurcation analysis.There exist numerous novel dynamical phenomena whose existence is fundamentally relatedto the presence of symmetry [16, 29, 31, 26]. Solutions related to the presence of symmetry,include rotating waves [54], modulated waves [49, 2], slow “phase” drifts along directionsof broken symmetry [41], and stable heteroclinic cycles [2, 33, 21]. The SO(2) symmetryis present in problems described by the Navier-Stokes [30, 62] or the Kuramoto-Sivashinsky[21, 3] equations with periodic boundary conditions, since the equations are invariant undertranslations and the boundary conditions do not break this invariance. Additionally, if thereflection symmetry exists, the full symmetry group is the O(2) group. While in classicalconvection problems (with constant viscosity), the study of the solutions and bifurcations inthe presence of symmetries has been the object of much attention [38, 32, 46, 43, 40, 39, 19, 4],its counterpart in fluids with viscosity depending on temperature has received less considera-tion. Our 2D physical set-up is idealized with respect to realistic geophysical flows occurringin the Earth’s interior, as these are 3D flows moving in spherical shells [11, 12]. Under theseconditions, the symmetry present in the problem is formed by all the orientation preservingrigid motions of R that fix the origin,which is the SO(3) group [14, 28, 36]. The effects ofthe Earth’s rotation are negligible in this respect and do not break this symmetry, since thehigh viscosity of the mantle makes the Coriolis number insignificant. The link between oursimplified problem and these realistic set-ups is that the O(2) symmetry is isomorphic to therotations along the azimuthal coordinate, which form a closed subgroup of SO(3). Addition-ally, the O(2) symmetry is present in systems with cylindrical geometry, which provide anidealized setting for volcanic conduits and magma chambers. The SO(2) symmetry is alsopresent in 3D flows moving in spherical shells which rotate around an axis.In this way, specific symmetry-related solutions found in our setting are expected tobe present in these other contexts. The interest of 2D numerical studies for representing3D time-dependent thermal convection with constant viscosity has been addressed in [56].The authors report that in turbulent regimes at high Rayleigh numbers, the flow structureand global quantities such as the Nusselt number and the Reynolds number show a similarbehaviour in 3D and 2D simulations as far as high values of the Prandtl number are concerned.2igure 1: Problem set-up. A 2D container of length L and depth d with periodic lateralboundary conditions. The bottom plate plate is at temperature T and is rigid, while at theupper plate the temperature is T ( T < T ) and free slip.In some sense, these results suggest that our simulations might be illustrative for the 3Dcase, since although they are far from a turbulent regime and do not correspond to the caseof constant viscosity, they have been performed according to the infinite Prandtl numberapproach. In this article we show that typical solutions of systems with symmetries, aspreviously reported in diverse contexts [2, 21, 66], are also present in mantle convection andmagma-related problems. We report the presence of traveling waves and limit cycles nearheteroclinic connections after a Hopf bifurcation. We do this by means of bifurcation analysistechniques an direct numerical simulations and of the full partial differential equations system.The article is organized as follows: In Section 2, we formulate the problem, providing thedescription of the physical set-up, the basic equations and boundary conditions. In Section 3we present the viscosity law under consideration and discuss several limits in which previouslystudied dependencies are recovered. Section 4 summarizes the numerical methods used tosketch an outlook of the solutions displayed by the system. Section 5 discusses the solutionsobtained for a broad parameter set. Finally Section 6 presents the conclusions. As shown in Figure 1 we consider a fluid layer, placed in a 2D container of length L ( x coordinate) and depth d ( z coordinate). The bottom plate is at temperature T and theupper plate is at T , where T = T − ∆ T and ∆ T is the vertical temperature difference,which is positive, i.e , T < T .The magnitudes involved in the equations governing the system are the velocity field u = ( u x , u z ), the temperature T , and the pressure P . The spatial coordinates are x and z andthe time is denoted by t . Equations are simplified by taking into account the Boussinesq ap-proximation, where the density ρ is considered as constant everywhere except in the externalforcing term, where a dependence on temperature is assumed, as follows ρ = ρ (1 − α ( T − T )) . Here ρ is the mean density at temperature T and α the thermal expansion coefficient.The equations are expressed with magnitudes in dimensionless form after rescaling asfollows: ( x (cid:48) , z (cid:48) ) = ( x, z ) /d , t (cid:48) = κt/d , u (cid:48) = d u /κ , P (cid:48) = d P/ ( ρ κν ) , θ (cid:48) = ( T − T ) / (∆ T ).Here, κ is the thermal diffusivity and ν is the maximum viscosity of the fluid, which isviscosity at temperature T . After rescaling the domain, Ω = [0 , L ) × [0 , d ] is transformedinto Ω = [0 , Γ) × [0 ,
1] where Γ =
L/d is the aspect ratio. The system evolves accordingto the momentum and the mass balance equations, as well a to the energy conservation3 ( ) / Ra=600Ra=800Ra=1000Ra=1200Ra=1600Ra=2000 (a) ( ) / b=1b=5b=10b=17b=23b=30 (b) ( ) / a=0.1a=0.01a=0.001 (c) Figure 2: Representation of the arctangent viscosity law versus the dimensionless temperaturefor different parameters values; a) b = 10, a = 0 . a = 0 . b values; c) b = 10, Ra = 1300 and different a values.principle. The non-dimensional equations are (after dropping the primes in the fields): ∇ · u = 0 , (1)1Pr ( ∂ t u + u · ∇ u ) = Ra θ(cid:126)e − ∇ P + div (cid:18) ν ( θ ) ν ( ∇ u + ( ∇ u ) T ) (cid:19) , (2) ∂ t θ + u · ∇ θ = ∆ θ. (3)Here, (cid:126)e represents the unitary vector in the vertical direction, Ra = d αg ∆ T / ( ν κ ) isthe Rayleigh number, g is the gravity acceleration, Pr = ν /κ is the Prandtl number. Typi-cally for rocks Pr, is very large, since they present low thermal conductivity (approximately10 − m /s ) and very large viscosity (of the order 10 N s/m ) [20]. Thus, for the problemunder consideration, Pr can be considered as infinite and the left-hand side term in (2) canbe made equal to zero. The viscosity ν ( θ ) is a smooth positive bounded function of θ , whichin our set-up represents a transition in the fluid, due for instance to the melting of mineralscaused by an abrupt change in viscosity at a certain temperature. This is discussed in detailin the following section.For the boundary conditions, we consider that the bottom plate is rigid and that theupper surface is non-deformable and free slip. The dimensionless boundary conditions areexpressed as, θ = 1 , u = (cid:126) , on z = 0 and θ = ∂ z u x = u z = 0 , on z = 1 . (4)Lateral boundary conditions are periodic. Jointly with equations (1)-(3), these conditions areinvariant under translations along the x -coordinate, which introduces the symmetry SO(2)into the problem. In convection problems with constant viscosity, the reflexion symmetry x →− x is also present insofar as the fields are conveniently transformed as follows ( θ, u x , u z , p ) → ( θ, − u x , u z , p ). In this case, the O(2) group expresses the full problem symmetry. The newterms introduced by the temperature dependent viscosity, in the current set-up equation (2)maintain the reflexion symmetry, and the symmetry group is O(2). We consider that the viscosity depends on temperature, and that it changes more or lessabruptly at a certain temperature interval centered at a temperature of transition. This isexpressed with an arctangent law which reads as follows: ν ( T ) = A arctan( β { ( T − T ) − b } ) + A (5)The parameter β controls how abrupt the transition of the viscosity with temperature is. Veryhigh β values imply that the viscosity transition occurs within a very narrow temperature4ap, while a finite and not too large value β assumes that the phase change happens over amushy region of finite thickness [64]. For the results reported in this article, we have fixed β = 0 .
9. As β is fixed, the viscosity transition always occurs in a temperature interval withconstant thickness ∆ θ ∼ .
23. The temperature at which the transition occurs is controledby b . The constants A and A are adjusted by imposing that at the reference temperature T the viscosity law (5) must be ν . On the other hand, in the limit T >> T , for instance T − T = 2500, the viscosity is fixed to a fraction a of the viscosity ν . These conditionssupply the system: ν = A arctan( − βb ) + A ν a = A arctan( β { − b } ) + A which has the solution: A = ν (1 − a )arctan( − βb ) − arctan( β (2500 − b )) ,A = ν − A arctan( − bβ ) . In dimensionless form, the viscosity law becomes: ν ( θ ) ν = C arctan( β (Ra θµ − b )) + C (6)where C = A /ν and C = A /ν . In this expression, Ra is the Rayleigh number, θ is thedimensionless temperature, which takes values between 0 at the upper surface and 1 at thebottom. The parameter µ , defined as µ = ν κ/ ( d αg ), is in this study fixed to µ = 0 . a is related to the inverse of the maximum viscosity contrast on the fluid layer,although the viscosity ν a may not correspond to any element of the fluid layer. For instanceFigure 2(a) shows the viscosity variation with temperature for different Rayleigh numbers at a = 0 . b = 10. It is observed that at low Ra, Ra = 600, the viscosity is almost uniformin the fluid layer, and it is only beyond Ra = 1000 that the sharp change in the viscosity isperceived. Figure 2(b) shows the effect of varying b at Ra = 1300 and a = 0 .
1. If b is as smallas 1, the transition occurs close to θ = 1 and most of the layer has low viscosity, while if b is very large at this Ra number most of the fluid has constant viscosity ν . It is interestingto relate the viscosity law as represented in these figures with the linear stability analysisof a fluid layer with constant viscosity ν , as presented in Figure 3. In this figure, one mayobserve that the critical Ra number is approximately Ra c ∼ b is large, the viscosity near the critical Rayleigh numberis almost constant across the fluid layer. In this case, the phase transition is noticed in thefluid at large Ra numbers, well above Ra = 1300, in a convection state in which vigorousplumes are already formed, as may be deduced from Figure 2(a). Figure 3(a) confirms thatat this limit the instability threshold of the conductive state remains unchanged with respectto that obtained with constant viscosity. On the other hand, if b is small, changes in the fluidviscosity are noticed at low Ra numbers –below the critical threshold of a fluid with constantviscosity– and in this case the instability threshold of the conductive state is affected bythe phase transition. This is illustrated, for instance, in Figures 2(a) and 3(b). For b = 10and a = 0 .
1, the changes in the viscosity across the fluid layer are noticed from Ra = 800onwards, which is below the instability threshold obtained for constant viscosity. In this case,the instability thresholds for the conductive solution are as those displayed in Figure 4, andthus the phase transition is perceived from the beginning by weakly convective states.We now discuss the relation between the arctangent law and an Arrhenius type lawfrequently used in the literature to model mantle convection problems. This viscosity law isexpressed according to [37, 20] as: ν ( θ ) = ν exp (cid:20) E ∗ R ∆ θ (cid:18) θ + t −
11 + t (cid:19)(cid:21) (7)5 R m=1 m=2 m=3 m=4 m=5 m=6 (a) R m=1 m=2 m=3 m=4 m=5 m=6 m=7 (b) Figure 3: Critical instability curves of the Rayleigh number, Ra, versus the aspect ratio Γat different wave numbers m . The results are for a fluid layer a) with constant viscosity; b)with temperature dependent viscosity µ = 0 . a = 0 . b = 30. ! R m=4m=1 m=7 m=8 m=9 m=11m=10m=6m=5m=3m=2 Figure 4: Critical instability curves of the Rayleigh number, Ra, versus the aspect ratio Γat different wave numbers m . The results are for a fluid layer with temperature dependentviscosity µ = 0 . b = 10 and a = 0 . a = 0 .
01 (thin line).6 ( ) / Arrhenius law Arctangent law b=1 Arctangent law b=5 Arctangent law b=10 Arctangent law b=30
Figure 5: The law of the viscosity dependent on the temperature used in [37] with viscositycontrast of factor 10 against the arctangent law (6) with parameters b = 1 , , ,
30, Ra =2500 and a = 0 . E ∗ is the activation energy, R is the universal gas constant, ∆ θ is the temperaturedrop across the fluid layer and t is the surface temperature divided by the temperature dropacross the layer. Figure 5 represents the viscosity (7) versus the dimensionless temperaturefor E ∗ R ∆ θ = 0 . t = 0 . b values are displayed. In this representation, one may observe the greatsimilitude between the Arrhenius law and the arctangent law for b = 1. At larger b values, thedecaying rate between viscosities is still similar to an Arrhenius law; however, temperatureintervals exist with approximately constant viscosities ν and ν a .One of the effects of the viscosity contrasts in the fluid motion is that if they are verylarge, as achieved for instance with the exponential or the Arrhenius law, they lead to astagnant lid convection regime [44, 58, 59], in which there exists a non mobile cap whereheat is dissipated mainly by conduction over a convecting flow. In [64, 17] a similar stagnantregime is obtained for a viscosity law similar to the one presented in this section. In oursetting, we have considered a free slip boundary at the top boundary, thus quiescence is notimposed. This condition enables us to consider spontaneous transitions from stagnant to nonstagnant regimes. Analysis of the solutions to the problem described by equations (1)-(3) and boundary con-ditions (4) is assisted by time dependent numerical simulations and bifurcation techniquessuch as branch continuation. As highlighted by [18, 48], the combination of both techniquesprovide a thorough insight into the solutions observed in the system. A full discussion onthe spectral numerical schemes used is given in [18]. For completeness, we now summarizethe essential elements of the numerical approach.
The simplest stationary solution to the problem described by equations (1)-(3) with boundaryconditions (4) is the conductive solution which satisfies u c = 0 and θ c = − z +1. This solutionis stable only for a range of vertical temperature gradients which are represented by smallenough Rayleigh numbers. Beyond the critical threshold Ra c , a convective motion settlesin and new structures are observed which may be either time dependent or stationary. Inthe latter case, the stationary equations, obtained by canceling the time derivatives in the7ystem (1)-(3) are satisfied by the bifurcating solutions. At the instability threshold ofthe conductive state, the growing solutions are periodic and correspond to sine or cosineeigenfunctions with wave number m . Figures 3 and 4 display the critical instability curvesfor different m values as a function of the aspect ratio. The new solutions depend on theexternal physical parameters, and new critical thresholds exist at which stability is lost,thereby giving rise to new bifurcated structures. These solutions are numerically obtainedby using an iterative Newton-Raphson method. This method starts with an approximatesolution at step s = 0, to which is added a small correction in tilde:( u s + ˜ u , θ s + ˜ θ, P s + ˜ P ) . (8)These expressions are introduced into the system (1)-(3), and after canceling the nonlinearterms in tilde, the following equations are obtained:0 = ∇ · ˜ u + ∇ · u s , (9)0 = − ∂ x ˜ P − ∂ x P s + 1 ν [ L ( θ s , u sx , u sz ) + L ( θ s )˜ u x + L ( θ s )˜ u z + L ( θ s , u sx , u sz )˜ θ ] , (10)0 = − ∂ z ˜ P − ∂ z P s + 1 ν [ L ( θ s , u sx , u sz ) + L ( θ s )˜ u x + L ( θ s )˜ u z + ( L ( θ s , u sx , u sz ) + Ra)˜ θ ] , (11)0 =˜ u · ∇ θ s + u s · ∇ ˜ θ + u s · ∇ θ s − ∆˜ θ − ∆ θ s . (12)Here, L ij ( i = 1 , j = 1 , , ,
4) are linear operators with non-constant coefficients, whichare defined as follows: L ( θ, u x , u z ) =2 ∂ θ ν ( θ ) ∂ x θ∂ x u x + ν ( θ )∆ u x + ∂ θ ν ( θ ) ∂ z θ ( ∂ x u z + ∂ z u x ) , (13) L ( θ ) =2 ∂ θ ν ( θ ) ∂ x θ∂ x + ν ( θ )∆ + ∂ θ ν ( θ ) ∂ z θ∂ x , (14) L ( θ ) = ∂ θ ν ( θ ) ∂ z θ∂ x , (15) L ( θ, u x , u z ) =2 ∂ θ ν ( θ ) ∂ x u x ∂ x + 2 ∂ θθ ν ( θ ) ∂ x θ∂ x u x + ∂ θ ν ( θ )∆ u x + ( ∂ x u z + ∂ z u x )( ∂ θ ν ( θ ) ∂ z + ∂ θθ ν ( θ ) ∂ z θ ) , (16) L ( θ, u x , u z ) =2 ∂ θ ν ( θ ) ∂ z θ∂ z u z + ν ( θ )∆ u z + ∂ θ ν ( θ ) ∂ x θ ( ∂ z u x + ∂ x u z ) , (17) L ( θ ) = ∂ θ ν ( θ ) ∂ x θ∂ z , (18) L ( θ, u x , u z ) =2 ∂ θ ν ( θ ) ∂ z θ∂ z + ν ( θ )∆ + ∂ θ ν ( θ ) ∂ x θ∂ z , (19) L ( θ, u x , u z ) =2 ∂ θ ν ( θ ) ∂ z u z ∂ z + 2 ∂ θθ ν ( θ ) ∂ z θ∂ z u z + ∂ θ ν ( θ )∆ u z + ( ∂ z u x + ∂ x u z )( ∂ θ ν ( θ ) ∂ x + ∂ θθ ν ( θ ) ∂ x θ ) . (20)The unknown fields ˜ u , ˜ P , ˜ θ are found by solving the linear system with the boundary con-ditions: ˜ θ = 0 , ˜ u = (cid:126) , on z = 0 and ˜ θ = ∂ z ˜ u x = ˜ u z = 0 , on z = 1 . (21)Then the new approximate solution s + 1 is set to u s +1 = u s + ˜ u , θ s +1 = θ s + ˜ θ, P s +1 = P s + ˜ P .
The whole procedure is repeated for s + 1 until a convergence criterion is fulfilled. In partic-ular, we consider that the l norm of the computed perturbation should be less than 10 − .The study of the stability of the stationary solutions under consideration is addressed bymeans of a linear stability analysis. Now perturbations are added to a general stationarysolution, labeled with superindex b : u ( x, z, t ) = u b ( x, z ) + ˜ u ( x, z )e λt , (22) θ ( x, z, t ) = θ b ( x, z ) + ˜ θ ( x, z )e λt , (23) P ( x, z, t ) = P b ( x, z ) + ˜ P ( x, z )e λt . (24)8he sign in the real part of the eigenvalue λ determines the stability of the solution: if it isnegative, the perturbation decays and the stationary solution is stable, while if it is positivethe perturbation grows over time and the conductive solution is unstable. The linearizedequations are:0 = ∇ · ˜ u (25)0 = − ∂ x ˜ P + 1 ν [ L ( θ b )˜ u x + L ( θ b )˜ u z + L ( θ b , u bx , u bz )˜ θ ] (26)0 = − ∂ z ˜ P + 1 ν [ L ( θ b )˜ u x + L ( θ b )˜ u z + ( L ( θ b , u bx , u bz ) + Ra)˜ θ ] (27)0 =˜ u · ∇ θ b + u b · ∇ ˜ θ + u b · ∇ θ b − ∆˜ θ + λ ˜ θ, (28)where the operators L ij are the same as those defined in equations (13)-(20). Equations (25)-(28) jointly with its boundary conditions (identical to (21)) define a generalized eigenvalueproblem.The unknown fields Y of the stationary (9)-(12) and eigenvalue problems (25)-(28) areapproached by means of a spectral method according to the expansion: Y ( x, z ) = (cid:100) L/ (cid:101) (cid:88) l =1 M − (cid:88) m =0 b Ylm T m ( z ) cos(( l − x ) + (cid:100) L/ (cid:101) (cid:88) l =2 M − (cid:88) m =0 c Ylm T m ( z ) sin(( l − x ) . (29)In this notation, (cid:100)·(cid:101) represents the nearest integer towards infinity. Here L is an odd numberas justified in [18]. 4 × L × M unknown coefficients exist which are determined by a collocationmethod in which equations and boundary conditions are imposed at the collocation points( x j , z i ) , Uniform grid: x j = ( j −
1) 2 πL , j = 1 , . . . , L ;Gauss–Lobatto: z i = cos (cid:18)(cid:18) i − M − − (cid:19) π (cid:19) , i = 1 , . . . , M ;according to the rules detailed in [18]. Expansion orders L and M are taken to ensureaccuracy on the results: details on their values are provided in the Results section. Together with boundary conditions (4), the governing equations (1)–(3) define a time-dependentproblem for which we propose a temporal scheme based on a spectral spatial discretizationanalogous to that proposed in the previous section. As before, expansion orders L and M are such that they ensure accuracy on the results and details on their values are given in thefollowing section. To integrate in time, we use a third order multistep scheme. In particular,we use a backward differentiation formula (BDF), adapted for use with a variable time step.The variable time step scheme controls the step size according to an estimated error E forthe fields. The error estimation E is based on the difference between the solution obtainedwith a third and a second order scheme. The result of an integration at time n +1 is acceptedif E is below a certain tolerance. Details on the step adjustment are found in [18].BDFs are a particular case of multistep formulas which are implicit , thus the BDF schemeimplies solving at each time step the problem (see [34]):0 = ∇ · u n +1 (30)0 = Ra θ n +1 (cid:126)e − ∇ P n +1 + div (cid:18) ν ( θ n +1 ) ν ( ∇ u n +1 + ( ∇ u n +1 ) T ) (cid:19) (31) ∂ t θ n +1 = − u n +1 · ∇ θ n +1 + ∆ θ n +1 , (32)where ∂ t θ n +1 is replaced by a backward differentiation formula.9n [18], it has been proved that instead of solving the fully implicit scheme (30)-(32),a semi-implicit scheme can produce results with a similar accuracy and fewer CPU timerequirements. The semi-implicit scheme approaches the nonlinear terms in equations (30)-(32) by assuming that the solution at time n + 1 is a small perturbation ˜ Z of the solutionat time n ; thus, z n +1 = z n + ˜ Z . Once linear equations for ˜ Z are derived, the equations arerewritten by replacing ˜ Z = z n +1 − z n . The solution is obtained at each step by solving theresulting linear equation for variables in time n + 1. In this section we explore how stationary solutions obtained at a low aspect ratio Γ = 3 . a and b of the viscosity law (6). We examinethe shape and structure of the plumes in a range of Rayleigh numbers from Ra = 2500 toRa = 3500.We first consider that the parameter b is large: for instance, as large as 30. In thiscase, Figure 2(b) confirms that at the instability threshold the viscosity across the fluidlayer is almost constant and equal to ν , no matter what the value of a may be. Thus,the viscosity transition becomes evident in the fluid once convection has settled in at Ranumbers well above the instability threshold. Figure 6(a) shows the plume pattern observedat Ra = 2500 for a = 0 .
1; although values a = 0 .
01 and a = 0 .
001 are not displayed,they provide a very similar output. The plume is spout-shaped, with the tail of the plumenearly as large as the head. In the pattern, the two black contour lines mark temperaturesbetween which the viscosity decays most rapidly. These correspond to the transition regionin which the gradient of the viscosity law (6) is large. Thus one of the contours, the coldestone, fits the temperature θ at which the viscosity has decayed by 5% from the maximum, i.e. , ν = 0 . ν , while the second addresses θ = θ + ∆ θ with temperature increment∆ θ = 0 .
23. The maximum viscosity decay rate always takes place at a constant temperatureincrement, since the decaying rate of the law (6) β , is the same through out all this study.At larger Rayleigh numbers, Ra = 3500, Figure 6(b) shows that the head of the plumebecomes more prominent. A comparison between Figure 6(b) and Figure 6(c) indicates thatthe large viscosity contrast favors the formation of a balloon-shaped plume, with a thinnertail and more prominent and rounded head. As regards the velocity fields, none of thesepatterns develop a stagnant lid at the surface for any of the viscosity contrasts a considered,even though the upper part corresponds to the region with maximum viscosity. This result isdissimilar to what is obtained in [64, 17]. In [17] it is argued that the cause of these differencescould be attributed to the transition sharpness controled by β , which in this work has beenconsidered to be smoother. Additionally, the results reported in [64] are obtained at largerviscosity contrasts, and the fact that these need to be large enough for the development of astagnant lid has been addressed.We now consider that the parameter b is small. As explained in Section 3, in this case theviscosity transition occurs at low Ra numbers, below the instability threshold of the fluid withconstant viscosity ν . As low viscosity also implies diminishing the critical Ra number, theoverall effect is that for small b the instability threshold is below that with constant viscosity ν , and the phase transition is perceived by weakly convective states. Figure 6(d) shows thestructure of the plume obtained for b = 10 and a = 0 . b = 5 or b = 1, exceptthat for smaller b values the tail of the plume tends to be thinner. Increasing the Ra numbermakes the tail of the plume thinner and spreads the head of the plume in the upper part, asreflected in Figure 6(e). On the other hand, high Ra numbers shift the viscosity transitiontowards colder temperature contours. As expected from the viscosity law (6), there is noRa number at which the whole fluid layer is “melted”, since this law always imposes that10 a) b = 30, Ra = 2500, a = 0 . b = 30, Ra = 3500, a = 0 . b = 30, Ra = 3500, a = 0 . b = 10, Ra = 2500, a = 0 . b = 10, Ra = 3500, a = 0 . b = 10, Ra = 2500, a = 0 . b = 17, Ra = 2500, a = 0 . b = 17, Ra = 2500, a = 0 .
001 (i) b = 17, Ra = 3000, a = 0 . Figure 6: (Color online). Plumes obtained for several values of the viscosity parameter b .The arrows indicate the velocity field, while the contour colors represent the temperatureranging from hot (bottom plate) to cold (upper plate). The two black contour lines indicatethe temperatures between which viscosity decays most rapidly.a transition occurs across the fluid layer. Figure 6(f) reports the effect of diminishing theviscosity contrast a to a = 0 .
001 at Ra = 2500. A mushroom-shaped plume with a thin tailand prominent head is observed. As before, none of these solutions develop a stagnant lid atthe surface for any of the examined viscosity contrasts a .Intermediate values such as b = 17 interpolate these extreme patterns. Figure 6(g) showsthe evolution from Figure 6(d) to Figure 6(a) in which the black contour lines indicating theposition of maximum viscosity decay converge towards the ascending plume boundary, thushighlighting its shape. The head of the plume shrinks and the tail strengthens. Diminishing a to the contrast 0.001 transforms the structure into a balloon-shaped plume (Figure 6(h)),while an increase in the Ra number spreads the head of the plume in the upper fluid towardsa mushroom-shaped plume.The structure of the observed plumes as a spout, balloon or mushroom shape followsthe schematic profiles reported in [37]. In the limit of low b , our viscosity law –as reportedin Section 3– converges towards the Arrehnius law used by these authors, and the plumeshapes reported there are similar to ours. However, a detailed comparison between bothworks is not possible as unlike these authors we include the Ra number in the viscosity law,since this provides a better expression of the realistic situation in which the increment of theRa number is performed by increasing the temperature differences between the bottom andupper surfaces. Other viscosity laws, such as the exponential law reported in [18] providedifferent plume structures, which are mainly spout shaped.The results reported in this section are obtained with expansions ( L × M = 37 × L × M = 47 × b = 10, a = 0 . .2 Bifurcation diagrams and time dependent solutions Solutions to the system (1)-(3) experience bifurcations depending on the aspect ratio andon the Rayleigh number. We now describe how these solutions vary along the dotted linesenhanced in Figure 4 for parameters µ = 0 . b = 10. We consider for a the choices0.1 and 0.01.Figure 7 shows the branch bifurcation diagram as a function of the aspect ratio forRayleigh number Ra = 1300 and a = 0 .
1. Branches are obtained by representing along thevertical axis the sum of the absolute value of two relevant coefficients in the expansion of thetemperature field, b θ and b θ . Solid lines stand for stable branches, while dashed lines arethe unstable ones. The horizontal line at | b θ | + | b θ | = 1 corresponds to the trivial conductivesolution. At a low aspect ratio, the stable branch is that with wave number m = 1, and ata higher aspect ratio the stable solutions increase their wave number to m = 2 and m = 3.The unstable branch ending up with a saddle-node bifurcation and connecting the m = 1with the m = 2 branch corresponds to a mixed mode.Stationary stable and unstable solutions, obtained at the positions indicated by arrows,are pictured. No stagnant lid appears at the surface for any of the aspect ratios considered.The expansion orders required by this figure to ensure accuracy are not the same along allbranches. We have guaranteed that for successive orders expansions the amplitude valuesdisplayed on the vertical axis of the bifurcation diagrams are preserved. A rule of thumb isthat high modes obtained at larger aspect ratios require higher expansions. Thus while formode m = 1 expansions ( L × M = 37 ×
44) are sufficient, for m = 2 and m = 3 at largeraspect ratios expansions are increased up to ( L × M = 61 × . a = 0 .
1. The pictured plumes, which are computed for a rather low Rayleigh number,Ra = 1500, are spout-shaped, with the tail of the plume nearly as large as the head. Asalready reported in the previous section for increasing Ra numbers, plumes become balloon-shaped and beyond that mushroom-shaped. No stagnant lid is observed at any Ra number.Several branches are distinguished. The branch related to mode m = 1 arises at the lowestRa number and is stable in the whole range displayed. Mode m = 2 emerges at Ra ∼ ∼ L × M = 37 × . a = 0 .
01. Figure9(a) examines the Ra interval from 800 to 1300. In this range several stationary solutionsare portrayed both stable and unstable. At Ra ∼ m = 3 (see Figure 9(b)). After the bifurcation, a traveling wave is found,as illustrated in the phase portrait represented at Ra = 1300. The solution evolves in timeby traveling towards the left. This breaks the symmetry x → − x . However, the righttraveling solution obtained by the symmetry transformation also exists, as expected fromequivariant bifurcation theory [16]. See [1] for further details. The presence of travelingwaves after a Hopf bifurcation has been reported in diverse contexts in under the presenceof the O(2) symmetry [2, 66, 21, 16], and here they are reported in the context of convectionwith variable viscosity. At larger Ra numbers, up to Ra ∼ m = 3 is found in therange Ra ∼ − ∼ ∼ ∼ b = 10, a = 0 .
1) at Γ = 3 . . Stationarysolutions are displayed at the Ra number, which is highlighted with the vertical line. Thearrows tag the branch points corresponding to the disclosed patterns.The dashed branchesare unstable, while the solid ones are stable. 14 a)(b)
Figure 9: (Color online). Bifurcation diagrams as a function of the Rayleigh number for afluid with viscosity dependent on temperature ( b = 10, a = 0 .
01) at Γ = 6 . . The dashedbranches correspond to stationary unstable solutions, while solid branches correspond tostationary stable ones. The gray lines indicate spatial patterns with period 2, while theblack ones are for period 3 patterns. a) Rayleigh number in the range 800-1300. Stationarysolutions are displayed at the Ra number highlighted with the vertical line. Arrows tagthe branch points corresponding to the disclosed patterns; b) Rayleigh number in the range1250-1500. Stationary solutions are displayed at the Ra number, which is highlighted bythe vertical line. The arrows tag the branch points corresponding to the disclosed patterns.Two additional vertical lines highlight the Ra = 1300 and Ra = 1400 numbers at whichtime dependent solutions are found. These are displayed as a time series projected on thecoefficient space (for a description see the text and [1]).15resence of plumes that are non-uniformly distributed along the horizontal coordinate: twoclose plumes, which are asymmetric around their vertical axis, and a third one that maintainsits symmetry. None of the described solutions develop stagnant lids at the surface. At low Ranumbers ( i.e.
Figure 9(a)) results are obtained with expansions ( L × M = 47 × , while forhigher Ra numbers ( i.e. Figure 9(b)) results are obtained with expansions ( L × M = 61 × . a = 0 .
1. The mode m = 3 branch, marked with a solid black line, emerges at Ra ∼ R ∼ R values marked with vertical dotted lines. A limit cycle is observed at Ra = 2210just above the bifurcation point. Its projection over the coefficient space displays a point atevery time step of the time series. The solution appears to reside in the neighbourhood of aheteroclinic connection between two fixed points as it evolves into a quasi-stationary regime–near the large density of points– followed by a rapid transition to a new quasi-stationaryregime. The two fixed points between which the solution oscillates are similar to the non-uniformly distributed plumes described in the previous paragraph (see [1] for further details).A solution is found at Ra = 2300 that has a time-dependence in which the block of plumesshifts irregularly along the horizontal direction, towards both the left and the right (see[1]). For increasing Ra numbers, the horizontal motion persists, but the oscillation becomesmore regular and pattern displacements along the x − coordinate are gradually reduced. Thisis verified through simulations at Ra = 2350 and at Ra = 2400 (see [1]). The diagramdisplayed in Figure 10(a) shows a gray solid line associated to a mode m = 2 stable branchthat emerges by means of a saddle node bifurcation jointly with an unstable branch. Anirregular pattern obtained at Ra = 1800 for the unstable branch is included in this diagram.Once again, none of the solutions described at this aspect ratio has a stagnant lid at thesurface. Results in this figure are obtained with different order expansions. At low R numberexpansions ( L × M = 47 ×
42) are sufficient while for higher Ra numbers they are increasedup to ( L × M = 61 ×
44) and even to ( L × M = 101 × This article addresses the study of a convection problem with temperature-dependent vis-cosity in the presence of the O(2) symmetry. In particular, the considered viscosity lawrepresents a viscosity transition at a certain temperature interval around a temperature oftransition. This is a problem of great interest for its many applications in geophysical andindustrial flows and in this work the focus is on exploring the impact of symmetry on thesolutions displayed by system.Our results report the influence on parameters a and b of the viscosity law on the mor-phology of the plumes at a low aspect ratio (Γ = 3 . ν , i.e , b is large, plumes tend to be thicker and show spout-like shapes. Increasing the Ra numberinduces their evolution towards balloon-shaped plumes, and this effect is more pronouncedfor high viscosity contrasts (small a ). At low b values plumes are thinner, and the head ofthe plume tends to spread in a mushroom-like shape in the upper part of the fluid.We explore bifurcations both for a fixed Ra number as a function of the aspect ratio, and16 a) | b | + | b | R Stables branches Unstable branches -0.1 0 0.10.010.0150.020.0250.03 b x c x -0.0500.05 -0.0200.02-0.1-0.0500.050.1 c x c x b x -0.0500.05 -0.0200.02-0.1-0.0500.050.1 c x c x b x -3 b x c x (b) Figure 10: (Color online). Bifurcation diagrams as a function of the Rayleigh number fora fluid with viscosity dependent on temperature ( b = 10, a = 0 .
1) at Γ = 7 . . The dashedbranches correspond to stationary unstable solutions, while solid branches correspond tostationary stable ones. The gray lines stand for spatial patterns with period 2 while theblack ones are for period 3 patterns. a) Rayleigh number in the range 700-1800. Stationarysolutions are displayed at the Ra number, which is highlighted by the vertical line. Thearrows tag the branch points corresponding to the disclosed patterns; b) Rayleigh numberin the range 1800-2500. Vertical lines highlight the Ra numbers at which time dependentsolutions are found. These are 2210, 2300, 2350 and 2400. These are displayed as a timeseries projected on the coefficients space (for a description see the text and [1]).17ifurcations at three fixed aspect ratios as a function of the Ra number. No stagnant lidregime is observed in any of the physical conditions analyzed. Among the stationary solutionsobtained along the bifurcation branches, one of the more interesting stable patterns consistsof the non-uniformly distributed plumes that break symmetry along their vertical axis.We also find that, for the higher Rayleigh numbers explored, at a high aspect ratio sev-eral rich dynamics appear. As already reported in classical convection problems, we finddynamical phenomena fundamentally related to the presence of symmetry, such as travel-ing waves, oscillating solutions in the neighborhood of heteroclinic connections and chaoticregimes characterized by “phase” drifts along the horizontal direction linked to the SO(2)symmetry.
Acknowledgements
We are grateful to CESGA and to CCC of Universidad Aut´onoma de Madrid for comput-ing facilities. This research is supported by the Spanish Ministry of Science under grantsMTM2008-03754, MTM2011-26696 and MINECO: ICMAT Severo Ochoa project SEV-2011-0087.
References [1] See supplemental material. http://link.aps.org/supplemental/10.1103/PhysRevE.88.043005 .[2] D. Armbruster, J. Guckenheimer, and P. Holmes. Heteroclinic cycles and modulatedtravelling waves in systems with o(2) symmetry.
Physica D , 29(257-282), 1988.[3] D. Armbruster, J. Guckenheimer, and P. Holmes. Kuramoto-sivashinsky dynamics onthe center-unstable manifold.
SIAM Journal on Applied Mathematics , 49(3):676–691,1989.[4] P. Assemat, A. Bergeon, and E. Knobloch. Nonlinear marangoni convection in circularand elliptical cylinders.
Physics of Fluids , 19:104101, 2007.[5] B. Blankenbach, F. H. Busse, U. Christensen, L. Cserepes, D. Gunkel, U. Hansen,H. Harder, G. Jarvis, M. Koch, G. Marquart, D. Moore, P. Olson, H.Schmeling, andM. Schnaubelt. A benchmark comparison for mantle convection codes.
Geophys. J. Int. ,98:23–38, 1989.[6] J. R. Booker. Thermal convection with strongly temperature-dependent viscosity.
J.Fluid Mech. , 76:741–754, 1976.[7] A. Bottaro, P. Metzener, and M. Matalon. Onset and two-dimensional patterns ofconvection with strongly temperature-dependent viscosity.
Physics of Fluids , 4(655-663), 1992.[8] G. Brandeis and C. Jaupert. On the interaction between convection and crystallizationin cooling magma chambers.
Earth Planet. Sci. Lett. , 77:345–361, 1986.[9] G. Brandeis and B. D. Marsh. The convective liquidus in a solidifying magma chamber:a fluid dynamic investigation.
Nature , 339:613–616, 1989.[10] H. P. Bunge, M. A. Richards, and J. R. Baumgardner. Effect of depth-dependent vis-cosity on the planform of mantle convection.
Nature , 436:436–438, 1996.[11] F. H. Busse. Pattern of convection in spherical shells.
J. Fluid Mech. , 72:65–85, 1975.[12] F. H. Busse and N. Riahi. Pattern of convection in spherical shells ii.
J. Fluid Mech. ,123:283–391, 1982. 1813] F. Capone and M. Gentile. Nonlinear stability analysis of convection for fluids withexponentially temperature-dependent viscosity.
Acta Mechanica , 107:53–64, 1994.[14] P. Chossat. Bifurcation and stability of convective flows in a rotating or not rotatingspherical shell.
SIAM Journal on Applied Mathematics , 37:624–647, 1975.[15] S. M. Copley, A. F. Giamel, S. M. Johnson, and M. F. Hornbecker. The origin of frecklesin unidirectionally solidified castings.
Metall. Trans. , 1:2193–2204, 1970.[16] J. D. Crawford and E. Knobloch. Symmetry and symmetry-breaking bifurcations influid dynamics.
Annu. Rev. Fluid Mech. , 23(341-387), 1991.[17] J. Curbelo and A. M. Mancho. Plate-like convection induced by symmetries in fluidswith temperature-dependent viscosity.
Preprint. arxiv: 1306.2921 , 2013.[18] J. Curbelo and A. M. Mancho. Spectral numerical schemes for time-dependent convec-tion with viscosity dependent on temperature.
Communications in Nonlinear Scienceand Numerical Simulations , 19, 2014.[19] P.C. Dauby, P. Colinet, and D. Johnson. Theoretical analysis of a dynamic thermocon-vective pattern in a circular container.
Phys. Rev. E , 61:2663, 2000.[20] G.F. Davies.
Dynamic Earth. Plates, Plumes and Mantle convection . Cambridge Uni-versity Press, 2001.[21] S. P. Dawson and A. M. Mancho. Collections of heteroclinic cycles in the Kuramoto-Sivashinsky equation.
Physica D: Nonlinear Phenomena , 100(3-4):231–256, 1997.[22] J. I. Diaz and B. Straughan. Global stability for convection when the viscosity has amaximum.
Continuum Mech. Thermodyn. , 16:347–352, 2004.[23] F. Dubuffet, D. A. Yuen, and M. Rabinowicz. Effects of a realistic mantle thermalconductivity on the patterns of 3d convection.
Earth Planet. Sci. Lett. , 171:401–409.,1999.[24] F. Dubuffet, D. A. Yuen, and T. K. B. Yanagawa. Feedback effects of variable thermalconductivity on the cold downwellings in high rayleigh number convection.
Geophys.Res. Lett. , 27:2981–2984, 2000.[25] J. D. Dufek and G. W. Bergantz. Transient two-dimensional dynamics in the upperconduit of a rhyolitic eruption: A comparison of the closure models for the granularstress.
J. Volcanol. Geotherm. Res , 143:113–132, 2005.[26] M. Field. Equivariant dynamical systems.
Trans. Am. Math. Soc. , 259(185-205), 1980.[27] G. R. Fulford and P. Broadbridge.
Industrial Mathematics . Australian MathematicalSociety Lecture Series 16. Cambridge University Press, 2002.[28] M. Golubitsky and D.G. Schaeffer. Bifurcation with o(3) symmetry including applica-tions to the benard problem.
Communs. Pure. Appl. Math. , 35:81–11, 1982.[29] M. Golubitsky and D.G. Schaeffer.
Singularities and Groups in Bifurcation Theory ,volume 1. Springer, 1985.[30] M. Golubitsky and I. Stewart. Symmetry and stability in taylor-couette flow.
SIAM J.Math. Anal. , 17(249-288), 1986.[31] M. Golubitsky, I. Stewart, and D.G. Schaeffer.
Singularities and Groups in BifurcationTheory , volume 2. Springer, 2nd edition, 2000.[32] M. Golubitsky, J. W. Swift, and E. Knobloch. Symmetries and pattern selection inrayleigh-benard convection.
Physica D , 10:249–276, 1984.1933] J. Guckenheimer and P. Holmes. Structurally stable heteroclinic cycles.
Math. Proc.Cambridge Philos. Soc. , 103(189-192), 1988.[34] E. Hairer, S.P. Norsett, and G. Wanner.
Solving Ordinary Differential Equations I.Nonstiff Problems . Springer, 2009.[35] J. M. Hyman, B. Nicoalaenko, and S. Zaleski. Order and complexity in the kuramoto-sivansinsky model of weakly turbulent interfaces.
Physica D , 23(1-3):265–292, 1986.[36] E. Ihrig and M. Golubitsky. Pattern selection with o(3) symmetry.
Physica D , 12:1–33,1984.[37] L. H. Kellogg and S. D. King. The effect of temperature dependent viscosity on thestructure of new plumes in the mantle: Results of a finite element model in a spherical,axisymmetric shell.
Earth and Planetary Science Letters , 148:13–26, 1997.[38] P. Kolodner, D. Bensimon, and C. M. Surko. Traveling-wave convection in an annulus.
Phys. Rev. Lett. , 60(1):723–726, 1988.[39] D. Krmpotic, B. Echebarria, and C. Perez-Garcia. Resonant interactions in benard-marangoni convection in cylindrical containers.
Physica D , 99(4):487–502, 1997.[40] D. Krmpotic, G. B. Mindlin, and C. Perez-Garcia. Benard-marangoni convection insquare containers.
Phys. Rev. E , 54(4):3609–3613, 1996.[41] M. Krupa. Bifurcations of relative equilibria.
SIAM J. Math. Anal. , 21:1453–1486, 1990.[42] L. G. Mastin. Insights into volcanic conduit flow from an open-source numerical model.
Geochem. Geophys. Geosyst. , 3:1037, 2002.[43] G. B. Mindlin, T. Ondarcuhu, H. L. Mancini, C. Perez-Garcia, and A. Garcimartin.Comparison of data from bernard-marangoni convection in a square container with amodel-based on symmetry arguments.
International Journal of Bifurcation and Chaos ,4(5):1121–1133, 1994.[44] L. N. Moresi and V. S. Solomatov. Numerical investigation of 2D convection withextremely large viscosity variations.
Physics of Fluids , 7(9):2154–2162, 1995.[45] M. Ogawa, G. Schubert, and A. Zebib. Numerical simulations of three-dimensionalconvection in a fluid with strongly temperature dependent viscosity.
J. Fluid Mech. ,233:299–328, 1991.[46] T. Ondarcuhu, G. B. Mindlin, H. L. Mancini, and C. Perez-Garcia. Dynamic patternsin bernard-marangoni convection in a square container dynamic patterns in bernard-marangoni convection in a square container dynamic patterns in bernard-marangoniconvection in a square container dynamic patterns in bernard-marangoni convection ina square container.
Phys. Rev. Lett. , 70(25):3892–3895, 1993.[47] E. Palm, T. Ellingsen, and B. Gjevik. On the occurrence of cellular motion in b´enardconvection.
J. Fluid Mech. , 30(651-661), 1967.[48] F. Pla, A. M. Mancho, and H. Herrero. Bifurcation phenomena in a convection prob-lem with temperature dependent viscosity at low aspect ratio.
Physica D: NonlinearPhenomena , 238(5):572–580, 2009.[49] D. Rand. Dynamics and symmetry: predictions for modulated waves in rotating fluids.
Arch. Ration. Mech. Anal. , 79(1):1–38, 1982.[50] J. T Ratcliff, P.J. Tacley, G. Schubert, and A. Zebib. Transitions in thermal convectionwith strongly variable viscosity.
Physics of the Earth and Planetary Interiors , 102:201–202, 1997. 2051] J. Revenaugh and B. Parsons. Dynamic topography and gravity anomalies for fluid layerswhose viscosity varies exponencially with depht.
Geophys. J. R Astrr. Soc. , 90:349–368,1987.[52] L. Richardson and B. Straughan. A nonlinear stability analysis for convection withtemperature–dependent viscosity.
Acta Mechanica , 97:41–49, 1993.[53] F. M. Richter, H. C. Nataf, and S.F. Daly. Heat transfer and horizontally averagedtemperature of convection with large viscosity variation.
J. Fluid Mech , 129:173–192,1983.[54] D. Ruelle. Bifurcations in the presence of a symmetry group.
Arch. Ration. Mech. Anal. ,51:136–152, 1973.[55] J. R. Sarazin and A. Hellawell. Channel formation in pb-sn, pb-sb and pb-sn-sb alloysand comparison with the system nh4cl-h2o.
Metall. Trans. , 19A:1861–1871, 1988.[56] J. Schmalzl, M. Breuer, and U. Hansen. On the validity of two-dimensional numericalapproaches to time-dependent thermal convection.
Europhysics Letters , 67(3):390–396,2004.[57] G. I. Sivashinky. On flame propagation under conditions of stoichiometry.
SIAM Journalon Applied Mathematics , 39(1):67–82, 1980.[58] V. S. Solomatov and L. N. Moresi. Stagnant lid convection on Venus.
J. Geophys. Res ,101:4737–4753, 1996.[59] V. S. Solomatov and L. N. Moresi. Three regimes of mantle convection with non-newtonian viscosity and stagnant lid convection on the terrestrial planets.
GeophysicalResearch Letters , 24(15):1907–1910, 1997.[60] V.S. Solomatov. Localized subcritical convective cells in temperature-dependent viscos-ity fluids.
Phys. Earth Planet. Inter. , 200-201:63–71, 2012.[61] V.S. Solomatov and A.C. Barr. Onset of convection in fluids with strongly temperature-dependent, power-law viscosity 2. dependence on the initial perturbation.
Phys. EarthPlanet. Inter. , 165(1-2):1–13, 2007.[62] I. Stewart and A. S. Hill. Three-mode interactions with o(2) symmetry and a model fortaylor-couette flow.
Dyn. Stab. Sys. , 6:267–339, 1991.[63] K. E. Torrance and D. L. Turcotte. Thermal convection wiith large viscosity variations.
J. Fluid Mech. , 47(113), 1971.[64] M. Ulvrov´a, S. Labrosse, N. Coltice, P. Raback, and P.J. Tackley. Numerical modelling ofconvection interacting with a melting and solidification front: Application to the thermalevolution of the basal magma ocean.
Physics of the Earth and Planetary Interiors , 206-207:51–66, 2012.[65] A. Vaidya and R. Wulandana. Non-linear stability for convection with quadratic tem-perature dependent viscosity.
Math. Meth. Appl. Sci. , 29:1555–1561, 2006.[66] S. A. van Gilsa and J. Mallet-Paret. Hopf bifurcation and symmetry: travelling andstanding waves on the circle.
Proceedings of the Royal Society of Edinburgh: Section A ,104(3-4):279–307, 1986.[67] D. B. White. The planforms and onset of convection with a temperature dependentviscosity.
J. Fluid. Mech. , 191:247–286, 1988.[68] T.K.B. Yanagawa, M. Nakada, and D.A. Yuen. A simplified mantle convection modelfor thermal conductivity stratification.