{'fraction_non_alphanumeric': 0.028845487914675533, 'fraction_numerical': 0.012016067594708775, 'mean_word_length': 4.764271457085829, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Symmetry and plate-like convection in fluids with temperature-dependent viscosity

We explore the instabilities developed in a fluid in which viscosity depends on temperature. In particular, we consider a dependency that models a very viscous (and thus rather rigid) lithosphere over a convecting mantle. To this end, we study a 2D convection problem in which viscosity depends on temperature by abruptly changing its value by a factor of 400 within a narrow temperature gap. We conduct a study which combines bifurcation analysis and time-dependent simulations. Solutions such as limit cycles are found that are fundamentally related to the presence of symmetry. Spontaneous plate-like behaviors that rapidly evolve towards a stagnant lid regime emerge sporadically through abrupt bursts during these cycles. The plate-like evolution alternates motions towards either the right or the left, thereby introducing temporary asymmetries on the convecting styles. Further time-dependent regimes with stagnant and plate-like lids are found and described.


Introduction
Rayleigh-Bénard convection is the classic example of thermal convection [2]. In these systems, under certain critical conditions, small fluctuations lead to massive reorganization of the convective motions [24,17,18,26]. This is a characteristic phenomenon of open systems that transfer energy and that are modelled by nonlinear equations. The internal energy of the planetary interiors is dissipated by convective processes, thus convection plays a crucial role in the evolution of the planet. Convective styles in planetary interiors are different from Rayleigh-Bénard convection. For instance, plate tectonics, which is distinctive of the Earth [21], is the surface manifestation of convection in the Earth's mantle. Other bodies in the solar system, such as the Moon, Venus or Mars do not exhibit plate tectonics and present other convection expressions with a stagnant lithosphere [32,33]. Different physical 1 justifications exist for the diverse types of convection: layered convection, for instance, is due to endothermic phase changes in the minerals that constitute the mantle interior [23]. The mantle is compressible due to changes in density, which increases towards the Earth interior. Numerical analysis of compressible convection indicates that density stratification has a stabilizing effect [20], produces upwelling plumes weaker than those downwelling and influences the thermal boundary layer [22] . The dependence of conductivity on temperature introduces new nonlinearities into the heat equation, which may lead to diverse dynamics [16]. When conductivity decreases with temperature, convection becomes more chaotic and time-dependent [12,37]. Thermal conductivity variation has generally been less studied than that of viscosity, as the latter is much stronger in the Earth's mantle. Large viscosity contrasts in fluids with temperature-dependent viscosity lead to stagnant lid convection [25,31]. Regarding the subduction initiation, numerical results [34,28,3,30] suggest that this is only possible if the stiff upper layers of the lithosphere are weakened by brittle fracture. Several mechanisms have been proposed for driving the motion of the lithospheric plates. Forsyth and Uyeda [13], for instance, conclude that plate-like motion is produced by the sinking slab that pulls the plate in the subduction process due to an excess of lithosphere density.
Finding the impact of the different physical properties present in the mantle on its convection styles is an important goal of research into planetary interiors. In this context, our focus is on examining the instabilities found in a 2D fluid in the presence of the O(2) symmetry which contemplates a phase transition similar to a melting-solidification processes. In particular, we consider a highly viscous layer (lithosphere) over a fluid mantle which is modeled with a viscosity that changes abruptly by a factor of 400, in a narrow temperature gap at which magma melts. In phase transitions, other fluid properties in addition to viscosity may change abruptly, such as density or thermal diffusivity. However, in this study we confine ourselves solely to the effects due to the variability of viscosity, since consideration of the effect of simultaneous variations on all the properties prevents a focused understanding of the exact role played by each one of these properties. Viscosity is a measure of fluid resistance to gradual deformation, and in this sense highly viscous fluids are more likely to behave rigidly when compared to less viscous fluids. When examining the proposed transition with temperature, we focus on the global fluid motion when some parts of this motion tend to be more rigid than others. By disregarding the variations on density in this transition, we move away from 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á et al. [35] deals with a problem similar to the one we address here, but takes into account variations in both density and viscosity. Thermal conductivity effects are related to the relative importance of heat advection versus diffusion. Diffusive effects are therefore important at large conductivity, while heat advection by fluid particles is dominant at low conductivity. The contrasts arising from these variations are beyond the scope of our work and are thus disregarded herein.
In our setting we show that convective processes exist which include plate-like motions that alternate in time with stagnant-lid regimes. Some of these transitions include bursts in which the solution releases energy to accommodate different spatial patterns. These solutions are mathematically related to limit cycles, which are persistent solutions in the presence of the O(2) symmetry [1,15,11] which is also found in this problem. There exist numerous novel dynamical phenomena in fluids that are fundamentally related to the presence of symmetries [7]: these include rotating waves [29], modulated waves [1,27] and stable heteroclinic cycles [1,15,11]. The SO(2) symmetry is present in the problem under consideration, because the equations are invariant under translations and periodic boundary conditions do not break this invariance. Additionally, if the reflection symmetry exists, the full group of symmetry is the O(2) group.
The impact of the symmetry on the solutions displayed in convection problems with temperature-dependent viscosity has been addressed in [9,8], where a 2D physical set-up similar to ours is analyzed. The viscosity law considered in this work is similar to the one studied in [8], the main difference being that the viscosity change in our current setting is achieved within a narrower temperature gap. Our problem is idealized in terms of realistic geophysical flows occurring in the Earth's interior, as these are 3D flows moving in spherical shells [4,5]. Under these conditions, the symmetry present in the problem is formed by all the orientation, preserving rigid motions of R 3 that fix the origin, which is the SO(3) group [6,14,19]. The effects of the Earth's rotation are negligible in this respect and do not break this symmetry, as the high viscosity of the mantle renders the Coriolis number insignificant. The link between our simplified problem and these realistic set-ups is that the O(2) symmetry is isomorphic to the rotations along the azimuthal coordinate, which form a closed subgroup of SO(3). Furthermore the O(2) symmetry is present in systems with cylindrical geometry, which provide an idealized setting for volcanic conduits and magma chambers. The results described in this paper confirm the symmetry role in the solutions that under the physical conditions considered exhibit plate-like dynamics and energy bursts.
The article is organized as follows. Section 2 describes the physical set-up and provides the governing equations as well as a detailed characterization of the viscosity law. Section 3 briefly introduces the numerical methods used to obtain the solutions. The results are presented in Section 4. Finally, Section 5 details the conclusions.

The physical set-up and the governing equations
We consider a convection problem in a fluid layer of thickness d placed in a 2D finite container of size L as shown in Fig 1. The bottom plate is rigid, i.e. u = 0, and it is at temperature T 1 . The upper plate is non-deformable and free slip and is at temperature T 0 , where T 0 = T 1 − ∆T and ∆T is the vertical temperature difference, which is positive, i.e, T 0 < T 1 . The lateral boundary conditions are periodic. The equations governing the system are expressed with magnitudes in dimensionless form after rescaling as follows: Here, κ is the thermal diffusivity, ρ 0 is the mean density at temperature T 0 and ν 0 is the reference viscosity. After rescaling the domain, 1] where Γ = L/d is the aspect ratio. The non-dimensional equations are (after dropping the primes in the fields): Here, e 3 represents the unitary vector in the vertical direction; Ra = d 4 αg∆T /(ν 0 κ) is the Rayleigh number; g is the gravity acceleration; α the thermal expansion coefficient and P r = ν 0 /κ is the Prandtl number. Typically for rocks, P r is very large, since they present low thermal conductivity (approximately 10 −6 m 2 /s) and very large viscosity (of the order 10 20 N s/m 2 ) [10]. Thus, for the problem under consideration, P r can be considered as infinite and the left-hand side term in (2) can be made equal to zero. These equations use the Boussinesq approximation in which the density is considered constant everywhere except in the buoyant term of Eq. (2) where a dependence on temperature is assumed, as follows ρ = ρ 0 (1 − α(T − T 0 )). Thus, no change in the density at the melting temperature is considered and is assumed to be small. Jointly with equations (1)-(3), the lateral periodic conditions are invariant under translations along the x-coordinate, which introduces the symmetry SO(2) into the problem. 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 symmetry of the problem.
The viscosity ν(θ) is a smooth, positive and bounded function of θ. We use a law that represents the melting by means of an abrupt change in the viscosity at a small temperature gap defining the melting transition. In dimensionless form, this law is, Here, the temperature at which the transition occurs is adjusted by the transition Rayleigh Ra t which in our case is Ra t = 10. The choice of a positive value for Ra t imposes that there exists a viscosity transition in the interior of the fluid layer, even if Ra is very large. The parameter β controls the abruptness of the viscosity transition on θ. Throughout this study we take β = 100. The constant µ, fixed to µ = 0.0146, expresses fluid properties. The presence of the Ra number in the viscosity law is uncommon among the literature dealing with viscosity dependent on temperature. However, it expresses better what happens in laboratory experiments in which the increment of the Ra number is performed by increasing the temperature T 1 at the lower boundary. This procedure ties the viscosity to changes in the Rayleigh number, which is the parameter that we vary in our study. Changes in the Ra number, as it appears in the viscosity law, necessarily imply changes in the viscosity contrasts. This is explored in further detail. The maximum viscosity in the fluid layer is ν 0 , and this is the viscosity value used to define the dimensionless Rayleigh number Ra in Eq. (2). In practice, ν 0 is a viscosity only taken by the fluid at the upper surface where θ = 0 (i.e, at temperature T 0 ), in the limit of large Ra t and small Ra. The parameter a is related to the inverse of the maximum viscosity contrast and is fixed at 10 −3 . In the range of Ra numbers considered in this study, we obtain viscosity contrasts of the order of 3 · 10 2 − 4 · 10 2 . Figure 2 represents the viscosity law at different Ra numbers within the range considered in this work. It is observed that as the Ra number increases the viscosity transition occurs in a temperature gap closer to θ = 0, and therefore closer to the upper surface. The conductive solution (i.e. u c = 0) to the problem described by equations (1)-(3) corresponds to the linear temperature θ c = −z + 1. Fig. 1 shows two viscosity profiles as a function of the depth z for this particular temperature solution. These profiles are obtained at the same Ra numbers as in Figure 2, i.e. Ra = 50, 150, and they confirm that the viscosity transition occurs close to the upper surface. A viscosity law similar to that expressed in Eq. (4) has been proposed in [35,8]. The study in [35] also considers a change in the density at the temperature of transition, while in [8] the viscosity change occurs within a broader temperature gap.

Numerical methods
The results presented in this work are obtained by solving the basic equations and boundary conditions with the numerical techniques reported in [9]. Our analysis is assisted by timedependent numerical simulations and bifurcation techniques such as branch continuation. These schemes are briefly described below.  The simplest stationary solution to the problem described by equations (1)-(3) and their boundary conditions is the conductive solution that satisfies u c = 0 and θ c = −z + 1. This solution is stable only for a range of vertical temperature gradients that are represented by small enough Rayleigh numbers. Beyond the critical threshold Ra c , a convective motion settles in and new structures are observed, which may be either timedependent or stationary. In the latter case, the stationary equations, obtained by cancelling the time derivatives in the system (1)-(3) are satisfied by the bifurcating solutions. At the instability threshold of the conductive state, the growing solutions are periodic and correspond to sine or cosine eigenfunctions with wave number m. Figure 3 displays the critical instability curves for different m values as a function of the aspect ratio. These curves are obtained by means of a simplified linear stability analysis for the conductive solution, as reported in [9]. At the instability threshold around Ra ∼ 55, the viscosity law indicates (see Figure 2) that the viscosity transition takes place at the lowest temperatures across the fluid layer, which is near the fluid surface. Thus, at this threshold the fluid consists of a highly viscous layer over a fluid that is not so viscous and is starting its convection.

Stationary solutions and their stability
Beyond the instability thresholds displayed in Figure 3, new branches of stationary solutions arise that evolve with the external physical parameters. There also exist new critical thresholds at which stability is lost, thereby giving rise to new bifurcated structures. These stationary solutions are numerically obtained by using an iterative Newton-Raphson method as reported in [9,8].
The study of the stability of the stationary solutions under consideration is addressed by means of a linear stability analysis. To this end, a field Y representing the unknown physical magnitudes is decomposed into its stationary solution Y b and a perturbationỹ as follows: The sign in the real part of the eigenvalue λ determines the stability of the solution: if it is negative, the perturbation decays and the stationary solution is stable, while if it is positive the perturbation grows over time and the stationary solution is unstable. For each unknown field expression, (5) is introduced into the system (1)-(3) and the equations are linearized inỹ, which are assumed to be small (see [9,8] for details). Together with their boundary conditions, the equations define a generalized eigenvalue problem. The unknown perturbation fieldsỹ of the linear equations are approached by means of a spectral method according to the expansion: In this notation, · represents the nearest integer towards infinity. Here, L is an odd number as justified in [9]. 4 × L × M unknown coefficients exist that are determined by a collocation method in which equations and boundary conditions are imposed at the collocation points, according to the rules detailed in [9]. Expansion orders L and M are taken to ensure accuracy on the results: details of their values are provided in the Section 4.

Time-dependent schemes
The governing equations (1)-(3) and their boundary conditions define a time-dependent problem for which we propose a temporal scheme based on a spectral spatial discretization analogous to that proposed in the previous section. As before, expansion orders L and M are such that they ensure accuracy on the results; details on their values are given in the following section. In order to integrate in time, we use a third order multistep scheme. In particular, we use a backward differentiation formula (BDF) that is adapted for use with a variable time step. Details on the step adjustment are found in [9]. BDFs are a particular case of multistep formulas which are implicit. In [9], it is reported that instead of solving the fully implicit scheme, a semi-implicit scheme is able to provide results with a similar accuracy and fewer CPU time requirements, and this is the method we employ to obtain the time-dependent solutions.

Results and discussion
Our study is focused on the solutions displayed by this system at a fixed aspect ratio Γ = 2.166. As the system is forced to transport more energy by increasing the Ra number, the conductive solution becomes unstable, and new convective solutions are observed. The bifurcation point for this primary event occurs at Ra ∼ 55. Furthermore, beyond this point a sequence of bifurcations occur when Ra increases, which is described below. Figures 4 and 5 show the bifurcating branches captured at different Ra ranges. In the diagrams, the horizontal axis represents the Ra number, while on the vertical axis the system state is represented by a scalar given by the sum of two coefficients in the expansion (6) of a stationary solution: This amplitude is related to the energy in the temperature field. The horizontal line at 0 corresponds to the conductive solution which is always a stationary solution of the system.    For instance, at Ra = 115 the solution for the solid grey branch is displayed, which corresponds to a pattern with one plume. Two white lines indicate the temperature contours at which the viscosity mostly decay. Figure 6 displays the horizontal component of the velocity versus the z-coordinate for this solution. The highlighted thin layer at the upper part confirms the existence of a stagnant lid at the surface. At Ra = 87 a two-plume solution in an unstable branch is represented. Several pitchfork bifurcations occur from which stable branches emerge. At Ra = 115 the pattern of a two-plume stable solution is shown; as in the one-plume solution, it has a stagnant upper surface. This stable branch undergoes a Hopf bifurcation at Ra ∼ 123, after which time-dependent solutions are found. A projection on the expansion coefficients space of this time-dependent solution is displayed at Ra ∼ 123. This solution consists of two plumes each slightly oscillating around their axis below a stagnant lid. After the Hopf bifurcation, the unstable branch undergoes a pitchfork bifurcation at Ra ∼ 130 and a stable branch emerges. This branch merges again with the unstable branch at Ra ∼ 140. The crossing of branches at Ra ∼ 135 is an effect of the projection taken and does not represent any transcritical bifurcation.
When the system is forced to transfer higher energy rates, further time-dependent solutions are observed. This is confirmed in Figure 5, which describes a bifurcation diagram similar to the previous one but at higher Rayleigh numbers Ra ∈ [143, 162]. The black branch bifurcates through a Hopf bifurcation at Ra ∼ 144 towards time-dependent regimes that coexist with the grey stable branch. Hopf bifurcations also occur for a different black stable branch at Ra ∼ 149.3, as well as for the grey solid branch at Ra ∼ 158. Figure 5 shows several stationary solutions: at Ra ∼ 148 it displays a two-plume structure obtained over the marked branch. This solution coexists with a time-dependent solution observed over the entire shaded region, of which we provide a projection of the time series obtained at Ra = 148, which is explained in detail below. An unstable stationary solution over this branch is displayed at Ra ∼ 155. The pattern is asymmetric, with the plumes more prominent outwards. At this Ra number, the structure of a stationary solution on the stable branch also exhibits asymmetry, but with the heads more prominent inwards. All the stationary solutions have a stagnant upper surface. We next describe the temporal evolution observed in the shaded region of Figure 5. A summary of what is obtained at Ra = 148 is given in Figure 7. The system evolves as a limit cycle in a type of evolution similar to that close to heteroclinic cycles, a typical object of systems with the O(2) symmetry [1,11]. Over time, the system stays close to two quiescent states, which are distinguished in Figure 7(a) by the zero upper velocity for long periods. These states are interrupted by bursts in which energy at the upper surface is abruptly released. One of the two quiescent states is represented in the center panel of Figure 7(b). In these almost stationary positions, the system presents two plumes that oscillate very slightly and have a stagnant lid at the upper surface. The two states are distinguished by the fact that the plumes are slightly shifted along the horizontal direction. Figure 7(c) represents the horizontal component of the velocity u x versus the depth z for a quiescent state at time t = 0.14482. In the center panel one may observe the stagnant upper lid. The quasi-stationary regimes are connected by rapidly evolving transitions which release the energy very rapidly. The solution during these crises has the interesting characteristic that in these episodes it consists of "plate-like" convective styles. By "plate-like" motion we refer to the fact that the stagnant lid at the upper surface drifts alternately towards the right or the left as a block. The first and third panels in Figure 7(c) show a moving upper layer by displaying the horizontal component of the velocity u x versus the depth z. A thin lid is observed which moves like a rigid body without internal shear either to the right or to the left. We have verified that variations of the velocity within this thin layer are below 0.05%, and for this reason it has the appearance of a moving plate. In these short time intervals (short when compared to the duration of the quiescent states), a meandering jet develops simultaneously below the drifting surface, in which sinking and upwards currents are observable (see Figure 7(b)).
It is clear that the lateral boundary conditions (i.e. the symmetry) are important for this behaviour because they allow upper drift motions. Moreover, the shifted inactive states between which the alternancy appears are possible only in this scenario. On the other hand, in a recent work by Ulvrová et al. [35] the authors study a law similar to ours, although symmetry effects are absent and this type of transition is not reported. Nevertheless, it should be noted that symmetry is not a sufficient condition for this kind of behavior: for instance, previous results discussed in [8], also in the presence of the O(2) symmetry, and in a similar setting to ours, make no reference to "plate-like" convection nor stagnant lids, although symmetries do exert an influence on the described solutions. The main difference between the setting considered in this work and that presented in [8]  14 is that the transition of the viscosity with temperature is less abrupt than in our study. Furthermore results reported in [9], also in the presence of the O(2) symmetry and viscosity according to exponential law, indicate the existence of a stagnant lid, although no impact of the symmetry on the time evolution is found. The bursting solutions obtained in this study are justified in the framework of the symmetry presence, but this does not imply that episodic solutions cannot be obtained in other settings in which the symmetry is absent. For instance, in the numerical work by [36], catastrophic events in Earth's mantle are reported. However, these bursting solutions are rather different to ours insofar as they do not connect almost quasi-stationary states and the solutions are not linked to any plate-like behaviour. The work by [21] is connected to episodic plate reorganizations, although the authors obtain time series in which transitions are not so dramatic as ours. However, in their problem they consider a square box with periodic boundary conditions, so they have the symmetry O(2)×O(2). Unfortunately, they provide no discussion on the impact of the underlying symmetry on their solutions.
Our simulations indicate that the solution described in Fig. 7 becomes non-attractive at Ra ∼ 149.5. At this point, an initial data starting near this regime evolves towards a periodic solution, a projection of which is shown in Figure 5. This motion consists of slightly asymmetrical plumes, each one rapidly vibrating around its central vertical axis. This is a time-dependent solution in which the surface fluid remains stagnant. At Ra ∼ 160, Figure  5 shows the projection on the coefficient space of a time series obtained for the asymptotic time-dependent regime of one plume. Figure 8 shows the time evolution in detail. The system evolves in a periodic motion in which the upper surface drifts alternately towards the right and the left, which is confirmed by a time series of the horizontal velocity component at the surface. Long quiescent states between the drift motions are no longer observed, only a continuously oscillating motion. In Figure 8(b), snapshots of the temperature and velocity fields obtained at times t = 0.14474, 0.14482, 0.145 show the plate-like motion. Figure 8(c) enforces this vision by displaying the horizontal component of the velocity u x versus the depth z at the same selected times. A thin upper lid which moves consistently is observed. As in the previously described plate motion, a meandering jet develops below the drifting surface in which sinking and upwards currents are observable (see Figure 8(b)). These time-dependent solutions are obtained with expansions L × M = 47 × 50.

Conclusions
In this paper we address the subject of a convecting fluid in which viscosity depends on temperature. We examine a dependency which models an abrupt change in the viscosity in a gap around a temperature of transition. We explore the space of solutions at a fixed aspect ratio by means of bifurcation diagrams and time-dependent numerical simulations. We find time-dependent convection in which the symmetry plays an important role. In particular, we describe limit cycles and time periodic solutions which are similar to others found in several contexts in the literature (see [1,15,11]) in the presence of the O(2) symmetry.
The time evolution during the limit cycles presents two peculiarities: first of all, they are bursting solutions that release energy abruptly in time and secondly plate-like convection is observed during the bursts. Additionally, time-periodic solutions are found that have a similar plate-like dynamic with a smoother time evolution. No plate-like dynamics have hitherto been observed in this type of convection problem. For viscosity dependencies according to the Arrhenius law, or its approach by means of an exponential law, no temporal transitions between stagnant lids and drifting lids have been reported. Recent studies by Ulvrová et al. [35], who use a law similar to ours, do not report this type of transitions either, although symmetry effects are not considered in their study.
Forsyth and Uyeda [13] propose that plate-like motion is produced by sinking slabs that pull the plates in the subduction process. The results reported in our study are obtained for constant density within the Boussinesq approximation, and provide convection examples of moving plates that coexist with subsurface upwards and downwards meandering jets, but without a proper subduction. Obviously, these examples do not rule out the existence of subduction in the Earth, but rather propose a role played by the symmetry which can be particularly illustrative for understanding convective styles of the Earth prior to subduction, or that of other planetary bodies.