On the competition of transverse and longitudinal modes of Marangoni convection in a three-dimensional layer of viscoelastic ﬂuid

Within the vast array of applications encompassed by viscoelastic ﬂuids, some lack of knowledge seems to aﬀect the non-linear behavior of Marangoni convection when its typical initial unicellular and steady state is taken over by more complex ﬂow conﬁgurations. These still hide a not-fully-understood competition of complex and diverse physical mechanisms that determine the prevailing macroscopic behavior. In the present study, relevant insights are sought from consideration of the classical diﬀerentially heated rectangular layer of liquid with adiabatic bottom and top free surface. It is shown that, for increasing values of the Marangoni number and/or the elasticity parameter, this problem oﬀers a multifaceted spectrum of diﬀerent outcomes depending on the non-trivial interplay established between two distinct categories of disturbances (transverse and longitudinal). These are studied using a diversity of model types in which some processes are on or oﬀ to discern selectively their eﬀect in the laminar state and their contribution to the evolution of the system towards chaos. The characteristic marks by which the ensuing elastic turbulence can be distinguished from the companion Kolmogorov counterpart are highlighted through analysis of the emerging scaling laws in the velocity spectrum and the sensitivity of these to the intensity of the driving force and the considered elasticity level. It is shown that these two forms of turbulence can coexist in the considered problem.


I. INTRODUCTION
Thermocapillarity can be defined as the tendency or ability of liquids to develop surface flow in the presence of temperature gradients. It naturally stems from the physical dependence of surface tension on temperature for many known substances 1 . This dependence typically leads to the existence of a surface stress imbalance, which results in the emergence of surface flow and, eventually, ensuing bulk fluid motion driven by viscosity 2, 3 .
These effects are omnipresent in nature and technology; inorganic and organic material solidification, crystal growth from the melt, soldering, film processing are just examples of the myriad technological processes where a liquid is exposed to a gas in the presence of temperature gradients with various orientations.
These systems have gained considerable momentum in the recent years owing to the rising tide of new studies aiming to deeply investigate them after filtering out the concurrent influence of gravitational effects (this being made possible by sounding rockets of various types and the advent of space platforms 4,5 such as the International Space Station). They have also been addressed extensively through numerical approaches of various kinds (see  just to cite some relatively recent contributions). Despite the remarkable efforts of research groups with various backgrounds and objectives, and the different directions from which the problem has been tackled, however, as a comparative review of the literature would confirm, the physical understanding of surface-tensiondriven convection is far less advanced than of buoyancya) Electronic mail: marcello.lappa@strath.ac.uk driven convection 9 .
This realization requires a short excursus into what is known, the "advances", and the new questions, both general and system-specific, that recent progresses in turn suggest. The simplest way to enter such discussion, perhaps, is to consider that different mechanisms can be enabled in surface-tension driven flows depending on the specific situation considered. A first set of influential factors is represented by the considered fluid Prandtl number (P r defined as ratio of the fluid kinematic viscosity and thermal diffusivity), the geometry of the physical domain hosting the liquid and the relative direction of the dominant temperature gradient with respect to the interface. The last aspect has largely been used in the literature to set a distinction between two main categories of flows, namely that of Marangoni-Bénard (MB) convection (in which the temperature gradient is perfectly perpendicular to the surface) and that of classical Marangoni convection, where the imposed temperature difference is in the same direction as the interface.
This apparently innocuous detail has proven to lead to remarkable differences in the mechanisms driving fluid flow and the related hierarchy of bifurcations along the path that leads to completely chaotic states. In the former case, fluid flow is typically produced from an initially quiescent and thermally diffusive state as soon as the temperature difference (or the related Marangoni number) exceeds a given threshold. The details of the effectively emerging convective states depend on the considered fluid. In infinite layers or cavities with large (lengthto-depth) aspect ratio filled with a high-Pr liquid (e.g. an oil), the pattern takes on a specific topology corresponding to the honeycomb symmetry, i.e. the cells look like perfect hexagons with fluid reaching the interface at the center and sinking at the periphery. This morphology is retained if liquid metals (P r << 1) are considered in place of oils. In this case, however, the relative direction of ascending and descending currents is inverted, i.e. the fluid sinks at the center and moves towards the interface along the sides of the hexagons. Moreover, extended horizontal rolls can also appear and coexist with the inverted hexagons. It is only by a further increase in the Marangoni number that these steady flows can be taken over by oscillatory states or patterns with a different structure 10,11 .
If the direction of the temperature gradient is rotated by 90 • , the scenario changes completely. Exceeding a threshold in order to induce convection from initially motionless conditions is no longer needed (pervasive fluid motion is created as soon as a temperature gradient is established, no matter how small). This typically leads to perfectly parallel flow (aligned with the imposed temperature gradient) in infinite layers or, similarly, a single extended convective roll with axis perpendicular to the imposed temperature difference in slender finite-size systems (cavities). On increasing the Marangoni number, these initially steady convective states can support the onset and propagation of travelling waves with orientation depending on the considered value of the Prandtl number 12 ; namely, the disturbances have a relatively small or significant spanwise (in a direction perpendicular to the temperature gradient) component according to whether P r > 1 (e.g. oils 13 ) or P r < 1 (liquid metals and semiconductor melts 14 ). In the former situation, the waves essentially behave as disturbances spreading in the upstream direction, i.e. they propagate from the cold to the hot side of the considered cavity, while the liquid located in proximity to the free interface moves in the opposite direction [15][16][17] .
Another influential factor is represented by the "elasticity" of the fluid, i.e. the ability to retain stresses even in the absence of velocity gradients. If present, this quality or property can enable a significant departure from the scenario known for standard Newtonian fluids. For MB convection, in particular, the typical "signature" of the elastic effects is the occurrence of "overstability", i.e. the ability of these systems to transition directly from a quiescent state to oscillatory convection without passing through the intermediate stage of steady flow 18 . This specific change in the observable behavior generally occurs in conjunction with significant lowering of the required value of the Marangoni number in comparison to the equivalent (in terms of viscosity and thermal diffusivity) Newtonian liquid. Moreover, the perfect symmetry of hexagonal cells seen in the Newtonian case is lost in favor of less regular (essentially more chaotic) patterning behavior and "multiple" solutions (i.e. states that exist in parallel in the space of parameters, and can selectively be accessed according to the considered initial conditions 19 ). Other works of relevance to the subject include those by Siddheshwar et al. 20 , Ramkissoon et al. 21 , Hernández Hernández and Dávalos-Orozco 22 , and the recent review by Lappa 23 , where buoyancy convection was also considered.
Studies of such a kind are just the beginning. Looking forward, viscoelasticity is currently being regarded as a spark for finding exotic states or new branches of solution or new mechanisms that, for now, remain undefined ideas in physicists' theories [24][25][26][27][28] .
Still fewer articles have been devoted to the case of Marangoni viscoelastic flow (temperature gradient parallel to the interface). Some guidance on this matter is provided by the investigations by Hu et al. [29][30][31] where it has been shown that for high-Pr fluids, elasticity can cause a swap in the direction of propagation of the aforementioned waves (enabled as a result of the primary Hopf bifurcation of the flow). More precisely, the disturbances start spreading in the downstream direction when the elasticity overcomes a given level, and eventually they are completely replaced by steady longitudinal modes of convection if a second threshold is exceeded (in such circumstances rolls with their axes parallel to the imposed temperature gradient are produced in place of the standard ones with perpendicular axes). More recently, Patne et al. 32,33 have extended this line of inquiry considering layers subjected to both parallel and perpendicular vertical temperature differences (inclined temperature gradient), thereby revealing a variety of stabilizing and destabilizing effects depending on the dominant mechanism.
These studies have opened up a new perspective on the subject by unvealing fundamental issues that are shared among the two areas of surface-tension-driven flows with vertical and parallel temperature gradients, which, therefore, can be said to unify their study. However, the surface of these specific problems has been barely scratched. The initially largely theoretical development of the subject, mostly based on the application of linear stability analysis (LSA) to ideally infinite systems, has not received yet a continuation with non-linear approaches. In particular, an open question is whether the modes of convection that have been specifically identified through LSA can be excited at the same time and eventually result in multiple solutions and ensuing peculiar convective states. The present study may be regarded as an attempt along these lines. In short, it follows and integrates the line of inquiry started by Lappa and Ferialdi 19 about the non-linear analysis of MB flow in viscoelastic fluids, by considering the equivalent problem in a layer of liquid with parallel temperature gradient.
The problem is tackled in the framework of direct numerical solution of the governing non-linear equations and results are provided for different levels of elasticity and distance from the critical conditions. In addition, relevant theoretical links are also explored with respect to recent literature (not necessarily related to surfacetension-driven flows) where the general interconnections between ordinary (inertial) and elastic turbulence have been considered. Schematization of the problem and locations where numerical probes have been placed to study the timedependent behavior of the flow

II. MATHEMATICAL MODEL
A finite 3D layer of viscoelastic fluid is considered. The layer is delimited at the sides and at the bottom by noslip walls, while the top boundary is modelled as an adiabatic and non-deformable free surface. Such a geometry can be characterized synthetically through its (dimensional) depth h and the aspect ratios A x = /h = 20 and A z = w/h = 10 for the x and z directions respectively. Obviously, the associated dimensional problem also requires that the fluid type and the applied temperature gradient (∇T ) are specified.
In the present work, the considered viscoelastic fluid is a mixture of a Newtonian solvent and a polymeric (or elastic material) solute having dynamic viscosity η s and η p respectively. The resulting fluid has therefore a total viscosity η 0 = η s + η p .
Accordingly, the dimensional balance equations for mass, momentum and energy in their non-linear and time-dependent form can be cast in compact form as: where t * is the time, u * is the velocity, T * is the temperature, p * is the pressure, ρ is the density of the fluid, andτ * is the extra-stress tensor accounting for the viscoelastic effects (the asterisk ( * ) is used to highlight the fact that those unknowns are in dimensional form; moreover, α is the fluid thermal diffusivity).
As this system of equations is typically solved to determine the three primitive variables u * , p * and T * , this automatically results in an important issue (which needs to be pinpointed suitably here), i.e., it implicitly makes evident that an additional equation is needed to linkτ * to the other unknowns (known as "constitutive relationship"). In practice, existing models essentially implement such a coupling betweenτ * and the fluid velocity u * (the reader being referred, e.g. to  for the Maxwell and Walter's liquid B paradigms, respectively). In the present work, in particular, we use the so-called Oldroyd-B constitutive model 41,42 .
From a historical point of view, the first solid theoretical underpinnings for such a paradigm emerged naturally out of the mathematics behind the theory of dumbbell molecules developed by Bird et al. 43 for polymer solution rheology. Since then, it has enjoyed a widespread use for the mathematical closure of the aforementioned set of balance equations. In particular, broad consensus exists in the literature that this model can adequately represent a wide spectrum of elastic solutions known as 'Boger fluids' (a relevant example being represented by a polyacrylamide solution in a maltose syrup/water mixture; the reader being also referred to the additional examples provided at the beginning of section IV). The hallmark of these liquids is their ability to exhibit an essentially constant viscosity over a wide range of shear rates 44 . Suffice it to say that for a variety of flows, this model has been successfully employed for the characterization of the novel instabilities that arise in viscoelastic fluids due to elasticity alone, or due to a combination of elastic and inertial effects (and which, therefore, have no counterparts in Newtonian fluids). Relevant examples include (but are not limited to) parallel plates flows, the lid-driven cavity, Taylor-Couette flows, Dean and Taylor-Dean Flows and cone-and-plate flows (the interested reader being referred, e.g., to Castillo Sánchez et al. 45 and Li and Khayat 44 for an exhaustive treatment).
Here we limit ourselves to mentioning that such lines of inquiry have demonstrated that, in the limit as the Reynolds number tends to zero (creeping flow regime), viscoelastic shear flows with curved streamlines undergo purely elastic instabilities and that these are essentially driven by normal stress differences. Although, the spatiotemporal properties of the emerging flow depend on the considered flow topology, the Oldroyd-B paradigm has proven to represent properly the physics of these purely elastic instabilities and led to a reasonable qualitative prediction of the conditions required for their onset. The Oldroyd-B is also at the basis of the various studies described to a certain extent in the introduction where surface-tension driven flows and related instabilities were considered 19,[29][30][31][32][33] . It has also been used for other types of thermal convection (see e.g. Refs 23, 46-48 and refer-ence therein).
From a practical standpoint, with the Oldroyd-B, the aforementioned link betweenτ * and the fluid velocity u * is formalized as a transport equation for the viscoelastic stress tensor, namely: where λ is the relaxation time, η p is the aforementioned dynamic viscosity of the polymer solute and (·) is the transpose operator of the quantity inside the parenthesis. The integration in time of this equation can therefore provide the required values ofτ * .
Given the problem being addressed in the present work, however, additional modeling is required to account for the presence of another category of stresses present in the fluid, i.e. the extra thermocapillary (Marangoni) stresses produced (at the interface separating the liquid layer from the external gaseous environment) as a result of temperature inhomogeneities.
Unlike the other balance equations, which can be regarded all as transport equations, the balance of these extra stresses reduces to a simple equality, namely where τ * = η s ∇u * + (∇u * ) +τ * is the total stress tensor, I is the unit tensor,n in the unit vector orthogonal to the free surface and σ T = −dσ/dT | T =T0 > 0 accounts for the decreasing behavior of the surface tension vs temperature in the linear relationship σ(T ) = σ 0 − σ T (T − T 0 ) 49-52 , where T 0 is the reference temperature and σ 0 = σ(T 0 ).
Although it introduces an important (vital) coupling between momentum and temperature (the latter would otherwise behave as a passively transported scalar quantity), this equation formally plays the role of boundary condition in the considered problem. It must therefore be considered together with the other relevant boundary conditions (namely no slip conditions along solid walls, constant temperature along the walls delimiting the layer along the x axis and adiabatic behavior for all the other boundaries).
Putting the problem in a "practical" shape finally requires some extra effort to replace the dimensional governing transport equations and the related boundary conditions with an equivalent set of non-dimensional (general) mathematical entities.
These can be obtained by re-scaling all the lengths with the depth of the layer h, the velocity with α/h, time with h 2 /α, pressure with ρα 2 /h 2 , temperature with ΔT = T h − T c and the viscoelastic stress with η 0 α/h 2 . Additional benefits stemming from this practice relate to the possibility to obtain relevant independent non-dimensional groups (characteristic numbers governing the considered problem) directly as coefficients, which appear in the equations as a result of the non-dimensionalization process itself: ϑ ∂τ ∂t + u · ∇τ +τ = ζξ ∇u + (∇u) These non-dimensional numbers are the Prandtl number P r = ν 0 /α (ν 0 = η 0 /ρ), the elasticity number ϑ = λα/h 2 , the solvent-to-total viscosity ratio ξ = η s /η 0 , the viscosity ratio ζ = η p /η s (= (1 − ξ)/ξ), and the Marangoni number Ma = σ T ΔT h/η 0 α.

III. NUMERICAL METHOD
Experimental struggles with understanding the basic mechanism of viscoelastic flows make the empirical works on these problems extremely difficult/case dependent, which provides an indirect justification for the widespread habit of using numerical simulations to get relevant insights. The present study should be regarded just as another effort along these lines. As suitably pointed out in the present section, however, the effective implementation of a relevant numerical method is not as straightforward as one would imagine.
Although, it is primarily the physical relevance of the Oldroyd-B model to real elastic liquids (i.e polymer solutions) which has made it a model of choice here, special care must be put on the implementation of relevant countermeasures to keep the time-marching algorithm stable and prevent it from developing unphysical solutions.
For the present study, in particular, the OpenFOAM computational platform has been used to solve the overall set of governing equations (presented in section II). More specifically, the equations have been discretised in space and integrated in time using a control volume method and a segregated numerical procedure, respectively. The latter is a well-known numerical realization of the PISO algorithm originally elaborated by Issa 53 (see also Refs 54 and 55), with a collocated disposition for the primitive variables and the Rhie and Chow 56 interpolation scheme used to avoid pressure checherboarding. The reader specifically interested in this approach and the related underlying rationale may consider Lappa and Boaro 46 . For the sake of brevity, in the present section we limit ourselves to describing the specific strategies put in place to avoid the divergence of the simulations while retaining physical consistency.
For what concerns the viscoelastic model implementation, we have relied on rheoTool 57 , a versatile computational OpenFOAM toolbox, known for its ability 24,58-60 to keep the numerical integration process of the viscoelastic stress transport equation more stable (thereby allowing exploration of a wider region of the space of parameters). The underlying method relies on the logconformation tensor approach 58,61 . The related "principles" can be briefly illustrated by simply recalling that the viscoelastic stress tensorτ can be expressed as where A is the conformation tensor and I is the unit tensor. In general, the conformation tensor is positive definite, however, in proximity to a critical point (a "singularity" due to the hyperbolic nature of the viscoelastic stress transport equation, see Refs 62-66 ) the positiveness of A can drop thereby leading to divergence of the numerical procedure. If instead of considering A, its natural logarithm is used, Θ = ln(A) can remain definite positive and, accordingly, the aforementioned singularityrelated problem can be strongly mitigated. Put simply, to solve eq. (4), rheoTool re-writes it in terms of Θ and then, using an exponential transformation and eq. (10), it calculates the value of the stress tensor to be used in the momentum equation (Ref 58 and references therein).
In addition, to ensure high quality of the results, we have used a second order accurate backward scheme to discretize the equations in time, a second order accurate central difference scheme for the spatial discretization of the diffusive terms and a third order Cubista scheme for the analogous treatment of the convective terms. To avoid nonphysical oscillations, the Cubista scheme is implemented though a deferred correction approach and the non scalar quantities are handled in a component-wise way 58 .
Last but not least, attention has also been paid to the intrinsic singularities of Marangoni flow in finite-size geometries. These can manifest regardless of the elastic or non elastic nature of the considered liquid (i.e. also in Newtonian fluids 67,68 ) as they are not associated to specific equations, rather stem from the inconsistency of the velocity boundary conditions used for the free interface and the isothermal walls in the "corners" of the cavity (which explains why they are generally referred to as "viscous singularities"). For the convenience of the reader who is not an expert in Marangoni flow, a heuristic interpretation for their origin can be briefly provided as follows. On the solid boundary where the no-slip condition is applied, the viscous stress is 0, while, on the free surface, in general, it is not zero as the Marangoni stress condition is effective there (right hand side of eq. (5)). Due to such unphysical imbalance, a singularity is produced in the computational cell located between the interface and the perpendicular wall. While for low Prandtl number fluids (P r << 1) the singularity can be implicitly bypassed by the finite precision associated with the local approximation used for the discretization of the derivatives, for high-Prandtl number fluids it typically manifests itself causing an unlimited growth of the velocity in this cell as the cell size is reduced. As a result, a need arises for the physical consistency to be recovered, and, from a mathematical point of view, this can be achieved if the Marangoni stress is forced to vanish at the solid walls while keeping a continuous behavior. This is typically obtained by using a "regularization function".
In the present work, in particular, to "regularize" the flow (Refs 67 and 68), eq. (5) has been modified as follows: where, following Refs 67 and 68, the regularization function R(x) reads: where 0 ≤ χ ≤ 1 accounts for the percentage of regularized surface. For all the simulations presented in section IV, χ = 0.02.
The function R modifies the stress on a small percentage of cells located along the free surface, bringing its value on the first cell of the series (that in contact with the solid wall) to ≈ 0. In this way it filters out the viscous singularity of Marangoni flow. This artifice also guarantees the aforementioned continuity of the stress along the free surface. The validity of this approach has already been shown in a number of studies appearing in the literature 69 .

A. Validation
In order to demonstrate the reliability of the numerical algorithm described in the preceding section, a two-stage validation hierarchy has been implemented. To test the validity of the regularization strategy, we have considered a classical benchmark, i.e. the Hopf bifurcation of standard Marangoni flow in a Newtonian fluid. In particular, the results obtained with the present approach have been assessed against those provided by a repeatedly validated  in-house code, conceived to investigate similar fluid dynamics phenomena 23,68-71 .
The results of such an exercise (reported in TABLE I) show a good agreement between the flow-field frequencies. As extensively discussed in section I, hydrothermal waves (HTWs) originating from Marangoni instabilities propagate upstream. In FIG. 2 we report the temporal evolution of the HTW for a two-dimensional (2D) case with A x = 20, P r = 15, Ma = 3 × 10 4 (to be directly compared with those reported by Lappa 70 for the same values of the governing parameters).
As a second stage of such modus operandi (to verify separately the coherence of the viscoelastic kernel), we have considered the Hopf bifurcation corresponding to the overstable behavior of Marangoni-Bénard convection in an infinite layer.
In particular for an Oldroyd-B liquid with P r = 200, ξ = 0.1 and ϑ = 0.2 the bifurcation occurs for a value of the Marangoni number Ma cr ≈ 61 if the Biot number at the free surface is 1 72 . Here the Biot number is defined in the classical way as Bi = hL/k, where h is the convective heat transfer coefficient at the free surface, k is the thermal conductivity of the fluid and L is the characteristic length of the problem (here it indicates the depth of the layer).
To determine the disturbance growth rate, we have monitored the amplitude of the velocity signal provided by a numerical probe located in the centre of a 2D fluid domain having aspect ratio 7.5 (with cyclic conditions as- sumed at the lateral boundaries in order to mimic a system infinitely extended along one direction). The critical Marangoni number for the onset of oscillatory convection has been obtained by extrapolating the amplitude reported in TABLE II to A = 0. The resulting critical value has been found to be Ma cr ≈ 59, which is in a good quantitative agreement with the data available in the work by Lebon et al. 72 .

B. Mesh refinement study
Selecting a relevant mesh is a non-trivial task stemming from the opposite needs to achieve gridindependent results and minimize the computational cost as much as possible. This is generally achieved with abinitio coarse-grained numerical simulations for a representative case, followed by a progressive increase in the used numerical resolution (typically achieved by using more computational points) until the sensitivity of the results on the grid drops below a given (pre-fixed) threshold. There are therefore two levels of difficulties intrinsically embedded in such an exercise, one being the intrinsic computational 'overheads' of the considered problem (which determine the effective computational cost), another one being the issue of discerning what can really be considered "representative" in terms of convergence.
For what concerns the first aspect, we have already illustrated some of the typical computational drawbacks of this category of flows in the preceding two sections, where we have discussed the relevant countermeasures (to be typically implemented together with a choice of the time integration step much smaller than that required in the equivalent Newtonian fluid situation). The second aspect (the selection of a relevant case for the mesh sensitivity assessment) brings problems of its own, these being the problem "dimensionality" (i.e. two dimensions or three dimensions?), the flow regime (steady, oscillatory or turbulent) and, last but not least, the quantity or process to be monitored to draw some relevant conclusions.
As in the present work (as it will be illustrated in detail in section IV), the flow displays a marked tendency to produce three-dimensional (3D) states that differ significantly from the equivalent 2D counterparts, a quite obvious choice should be a fully 3D case. Since, as it will be yet illustrated in detail in section IV, the flow is also able to produce turbulent states for relatively small values of the Marangoni and elasticity numbers (and these states constitute integral part of our study), a representative case should therefore contain these elements too, i.e. be chaotic. This specific aspect, in turn, connects to the last of the issues mentioned above, i.e. the identification of relevant local or global parameters or trends to be checked as the numerical resolution is increased.
In order to satisfy the requirement about a flow that is at the same time 3D and chaotic, we have selected the case P r = 7, Ma = 500 and ϑ = 0.2. Accordingly (with the objective in mind of making the properties of the emerging solution somehow "quantifiable"), we have determined for each solution the related velocity power spectral density (PSD).
In particular, the spectra have been calculated for two different velocity (u y ) signals, one located in the center of the layer (probe 11 in FIG. 1) and one closer to the solid boundary (probe 12). FIG. 3 shows such spectra for four different meshes. The related legend displays the used number of points along x, y and z, respectively. In order to assess the overall level of convergence, we have also calculated a "global measure", i.e. the integral I of the PSD in the sub-range that follows the elastic turbulence decay law. This quantity is related to the total power carried in such a sub-range. These values are summarized in TABLE III.
As the Reader will realize by inspecting these data, moving from a mesh having 300 × 30 × 150 elements to one having 350 × 30 × 175 has a negligible effect on both the scaling law (we will be more precise about the im- plications of this observation later), the PSD amplitude, and the overall power carried in the spectra. For these reasons,a mesh with 300 × 30 × 150 elements has been adopted for the present analysis (able to guarantee reasonable grid-independence while limiting the otherwise excessive computational time).

IV. RESULTS
All the cases presented hereafter relate to an Oldroyd-B fluid having P r = 7 and ξ = 0.1 while the Marangoni (Ma) and the elasticity (ϑ) numbers are allowed to span relatively wide ranges ( Ma = 500, ϑ ∈ [0.03, 0.4] and Ma = 2500, ϑ ∈ [0.01, 0.065]). As already explained before, by using this constitutive model, we implicitly refer to a specific category of liquids known as "Boger fluids" 44 , relevant examples being represented by waterbased polymer dilute solutions with a relatively small amount of a polymer such a PAM, PEG, PEO, PVP or Xanthan Gum. Notably, at relatively high temperatures, e.g., water between 50 • C and 80 • C, the Prandtl number of such mixtures would be similar to that considered in the present work (which is also the value used by Li and Khayat 44 ). As studying the progression of the system from initial laminar conditions towards chaos is one of the main objectives here, the interval for the ϑ parameter for each value of Ma has been defined accordingly, i.e. we have increased ϑ until a turbulent state has been produced. For larger Ma this has been achieved for a smaller ϑ, which explains why the two selected ranges reported above for Ma = 500 and Ma = 2500 are different.
Prior to expanding on the results, we wish also to highlight that, following a common practice in the study of Marangoni flows (as also witnessed by the existing literature on this subject in the case of Newtonian fluids), the free surface of the layer has been considered adiabatic. This assumption has often been used in past studies to mimic the situation where the two side walls are heated and cooled 'symmetrically' with respect to the ambient temperature, i.e. their temperatures are T a + ΔT/2 and T a − ΔT/2, respectively, where T a is the ambient temperature. In such circumstances, the temperature of the free surface (it is almost uniform with the exception of the changes that occur in proximity to the walls) is almost identical to that of the gas ambient (thereby minimizing the interfacial heat exchange). As a next step of the mod-eling hierarchy undertaken in the present work, future efforts shall obviously be devoted to assessing whether the heat exchange due to asymmetrical heating can exert a stabilizing or destabilizing influence on such dynamics (as an example, for the case of Newtonian fluids, it is known that heat loss can induce complex, often non-monotonic, behaviors in the flow instability threshold 73,74 ).
Similarly, not to make the objective of the present study (and the corresponding space of parameters to be investigated) too wide or complex, as mentioned at the beginning of section II, the liquid layer aspect ratio is kept fixed (the investigation of this specific aspect, see, e.g., Refs 75,76 for the case of Newtonian fluids, being delayed to a future study).
Given this premise, the determination of the other dimensional quantities corresponding to the abovementioned non-dimensional parameters simply requires fixing relevant values for λ and α. By assuming λ = 10 −3 s and α = 1.4 × 10 −7 m 2 /s (typical realistic values for small polymer concentrations), in practice, the values given above for ϑ and Ma can be used to obtain the required depth of the layer and temperature difference, respectively. By means of this approach, it is easy to verify that the above range of ϑ related to Ma = 500 could be obtained with liquid layers having depth varying between 0.07 mm and 0.019 mm. Assuming as typical value for the derivative of the surface tension σ T = 10 −4 Nm −1 K −1 , the related temperature difference would span the interval from 10.3 • C to 37.6 • C.
In an analogous way, taking into account that the set of ϑ values considered here for Ma = 2500 would produce a maximum and a minimum depth of the liquid layer equal to 0.18 mm and 0.046 mm, respectively, the related temperature difference would vary between 29.6 • C and 76 • C.
Given the well-known ability of Marangoni flow to produce transverse or longitudinal modes of convection in Newtonian fluids depending on the value of the Prandtl number [12][13][14] , as well as in viscoelastic fluids 29-31 , a precise analysis hierarchy is implemented in the remainder of this work, by which these two categories of disturbances are partially separated.
In this regard, we wish to recall that, unlike the LSA approach, where the spatial structure of disturbances can be fixed "a priori", this is not possible when an approach like that implemented in the present study is used. Indeed, when the problem is solved in its nonlinear form, the disturbances are naturally produced and selected out of the full possible spectrum of allowed perturbations according to their growth rate, i.e. on a "mostdangerous disturbance" basis. A possible way to concentrate on a specific category of disturbances, therefore, simply consists of preventing the system from developing certain perturbations by forcing it to retain certain symmetries in space, e.g. by obliging it to maintain a two-dimensional behavior (2D simulations).
This approach is intentionally used here to distinguish 2D (transverse) disturbances from the more complex spectrum potentially excited when the flow is allowed to develop along all the spatial dimensions (3D simulations).

A. 2D disturbances
Although, as we will clearly show in section IV B, the modes of convection that spontaneously arise in viscoelastic thermocapillary convection are essentially 3D, in line with the theoretical predictions of available LSA studies [29][30][31] , disturbances that are essentially 2D in nature can also be expected for high-Pr fluids if Ma and ϑ are located in a certain subregion of the space of parameters. As outlined before, these can be captured and studied in great detail if the problem is investigated numerically under the constraint of two-dimensionality. Obviously, with 2D simulations, 3D perturbations are naturally filtered out thereby leaving the ground to transverse waves or other effects which develop along the direction of the interface.
Given such premises, it is also worth recalling (see again the arguments provided in the introduction) that, in general, no threshold in terms Ma must be exceeded in order to produce fluid motion. This is a typical feature of thermocapillary convection. For Newtonian fluids, if a temperature gradient is applied along the free surface, no matter how small, a single stationary circulation system is established inside the cavity in the form of a horizontally elongated roll. This is still true in a viscoelastic fluid provided both the Marangoni number and the elasticity number are sufficiently small. In particular, according to the present results, for the Marangoni number in the range 500 ≤ Ma ≤ 2500, a steady state is attained only if the elasticity number is smaller than a given ϑ cr (which depends on the considered value of the Marangoni number), whereas more complex phenomena are enabled as soon ϑ is increased beyond this limit.
Following a logical approach where situations of increasing complexity are presented as the discussion progresses, most conveniently, we start from the description of the results obtained for the branch of solutions corresponding to Ma = 500. In this case, according to the numerical simulations the value of the threshold corresponds to ϑ cr = 0.15 (as a low frequency disturbance arises and a wave becomes visible for this value of the elasticity parameter). This may be regarded as a metastable condition because the wave survives only for a limited period of time (after this time, an asymptotic state with the classical elongated single cell is recovered). For slightly larger values of the elasticity number, however, this disturbance is turned into a permanent wave traveling in the downstream direction and manifesting as a weak corrugation in the shape of the streamlines (FIG. 4c). This perturbation has clearly an elastic nature as witnessed by its connection to the elasticity number and its direction of propagation (as opposed to the classical hydrothermal wave in Marangoni flow, which, as shown in section III A) propagates upstream).
Notably, and in a quite unexpected way, a further increase in ϑ leads to a completely different situation. As the elasticity number exceeds a second critical threshold (ϑ cr2 ≈ 0.23), the original single-roll circulation system is taken over by a multicellular structure consisting of a train of cells, which spread continuously in the downstream direction. On further increasing this parameter, the behavior becomes chaotic as witnessed by the tendency of cells to merge randomly each other or, vice versa, split into smaller eddies (see FIG. 4d at t = 46.88 for ϑ = 0.4).
This remarkable change in the patterning behavior is reflected by an analogous modification visible in the velocity PSD (power spectral density P) distribution obtained for a signal probed in the centre of the layer (see FIG. 5). As qualitatively and quantitatively substantiated by this plot, the spectrum displays the typical properties of a turbulent flow, namely, a broad continuous power-law-decay region spanning more than two orders of magnitude in terms of frequency.
Another key observation concerns the tide displayed by such a distribution of amplitudes for ω > 10. Notably, it does not follow the classical Kolmogorov scaling law, i.e. P(ω) ∝ ω −5/3 , where ω is the angular frequency. Rather the exponent is ≈ −3.2 in the intermediate range of frequencies and −6.4 in the high-frequency interval, which indicates that the considered chaotic state is not inertial in nature (we will return to this very important concept later).
A larger Marangoni number (Ma = 2500) increases the range of values of the elasticity number in which a stable elastic wave can be obtained. The elastic disturbance still travels downstream, however, it can be obtained for a value of ϑ as low as 3.6 × 10 −2 , basically one order of magnitude smaller than that found for Ma = 500, which, in line with the information reported in section II about the genesis of elastic waves, we ascribe to the larger curvature of the streamlines of Marangoni flow for higher Ma.
As soon as ϑ exceeds the second threshold ϑ cr2 ≈ 4.5 × 10 −2 , distinct rolls become visible also for this Marangoni number.
The different stages of evolution taken by the flow in one period of oscillation for this value of ϑ can be seen in FIG. 6. Comparing this figure with the chaotic state reported in FIG. 4d, it becomes evident that the morphology and the number of the rolls is different. The angular frequency is ω = 91.6 and, interestingly, a small increase in the elasticity number has the effect of causing a decrease in such a value. Moreover, at ϑ = 4.575×10 −2 new frequencies start to populate the spectrum, thereby making the identification of a predominant frequency in the flow relatively difficult or meaningless.
The significance of the next two figures of the sequence (FIG. 7 and Figure 8) resides in their ability to reveal the progressive transitions undergone by the signals and the related spectrum, respectively, as the elasticity parameter is increased.
For ϑ = 4.5 × 10 −2 the signal exhibits a relatively regular behavior, and this is reflected by its spectrum (plotted on a logarithmic plane in FIG. 8a) where the aforementioned dominant frequency can be clearly recognized (ω = 91.6 together with its harmonics, i.e. higher order multiples). A small increase in ϑ leads to a more involved flow. Indeed, for ϑ = 4.575×10 −2 the signal displays in some intervals a regular and periodic behavior, alternated to small periods of chaotic busts. Finally, for ϑ = 6.5×10 −2 , a completely turbulent signal is produced and the same scaling already identified for Ma = 500 and ϑ = 0.4 can be recognized, namely, the "−3.2/−6.4" law.
To complement this scenario with additional relevant data, we have reported in FIG. 9 the frequency corresponding to the highest amplitude in the spectrum as a function of ϑ. The graph is terminated at ϑ 0.064 as beyond this point distinguishing clearly such a frequency is no longer possible. Such a plot is useful as the non-monotonic behavior evident there indicates that disturbances with different wavenumber (and therefore frequency) can be excited for the considered value of the Marangoni number and therefore compete to become the most energetic mode of convection in the flow. There is a continuous switch from one dominant mode to another and back to the original mode as the elasticity parameter is increased, which indicates that a discrete set of modes exist which are in competition. In this regard an analogy might be drawn with other viscoelastic phenomena where multiple solutions are known to be the rule rather than an exception 19,46 .
Interestingly, in FIG. 9 the curves seem to tend to an asymptotic state with ω ≈ 52, which indicates that the energy tends to reside on a smaller frequency as a completely chaotic state is approached. An explanation/justification for this trend can be elaborated in its simplest form on the basis of the two-fold argument that 1) the total energy of the system is fixed (let us recall that all these transitions are produced as a result of an increase in ϑ, not Ma, which is fixed to 2500) and 2) the system can take advantage of a larger set of scales to distribute its energy as the number of coexisting modes is increased (rather than requiring an increase in the frequency of the disturbances as one should expect for transitions produced by an increase in the Marangoni number).

B. 3D disturbances
After considering the problem in the simplified framework represented by a two-dimensional (2D) configuration, we now turn to interpreting the equivalent dynamics emerging when such a constraint is removed. This specific modus operandi obeys the logic illustrated at the beginning of section IV, i.e. the precise intention to discern intrinsically 3D effects through critical comparison of the results provided by 3D simulations with the "equivalent" 2D ones (i.e. conducted for the same set of parameters). As we will illustrate in detail in this section, a strategy based on the progression from 2D to 3D simulations is indeed instrumental in revealing the role played by the dimensionality of the problem (that is the number of space dimensions involved) and the different level of "criticality" of disturbances that break or do not break certain symmetries of the considered system.
For consistency with the 2D results presented earlier, such analysis is implemented starting from the lower end of the considered interval of Ma, namely, Ma = 500. The main outcomes of the related 3D explorative simulations are shown in FIG. 10, FIG. 12 and FIG. 11 for values of the elasticity number which decrease from a ϑ larger than that needed to produce time-dependence in 2D (ϑ cr2D = 0.15) to a value as small as 0.03 (namely, ϑ = 0.2, 0.15, 0.05, and 0.03 are considered).
The most striking outcome of these simulations is the evidence they provide about the ability of the effective problem dimensionality to cause a departure from idealized solutions (such as the elongated single transverse roll typical of stable Marangoni flow) even if relatively small values of ϑ are considered.
As a fleeting glimpse into FIG. 10 for ϑ = 0.03 would immediately confirm, if the processes that depend on the details of the third direction are not excluded, longitudinal rolls are formed, i.e. a new class of disturbances is enabled. Besides their different spatial structure (the axes of the rolls being essentially aligned with the x direction), these modes of convection also differ significantly from the waves typical of 2D flow due to their temporal behavior. A good impression of the overall unsteady three-dimensional motion associated with this kind of disturbances can be inferred from the different snapshots reported in this figure (FIG. 10) at different times.
Four initial parallel rolls can be spotted in FIG. 10a. These rolls maintain their position until a variation in their topology is produced in proximity to the heated boundary (FIG. 10b). Warmer fluid starts to rise inside the region initially occupied by colder and descending liquid. As a result of this effect, the colder region is split in two parts, thereby giving rise to a fork-like shape in the roll distribution (with a central jet of hot fluid surrounded by two parallel 'strips' of cold fluid, (FIG. 10c). Interestingly, this variation in the topological configuration of the rolls is not stationary, rather, in an initial stage this "defect" travels in the downstream direction (i.e. from the hot side towards the cold one, from left to right in the figure), thus taking a behavior formally similar to that observed for the elastic waves in the 2D case. The dislocation travels along the whole length of the cavity until it reaches the cold boundary (FIG. 10d). Here, as a result of its interaction with this wall, a new fork-like localized feature is created that is mirror symmetric with respect to the original one (this time a central jet of cold fluid surrounded by two parallel strips of hot fluid can be distinguished, FIG. 10e). At this stage, the feature carrying the colder fluid starts to travel in the upstream direction (thereby formally resembling the behavior of a classical hydro-thermal wave , FIG. 10f).
These findings are naturally complemented by those reported in FIG. 11a where the spatio-temporal map of the vertical component of the velocity (u y ) along a fixed line in the spanwise direction passing through the centre of the layer is shown (i.e. x = 10, y = 0.5, and 0 ≤ z ≤ 10). This map is particularly useful as it allows the derivation of additional properties of the pattern which would otherwise be hidden or less evident 47,77 . In particular, it reveals that the propagation in the upstream motion is faster than the downstream counterpart. The phenomenon repeats itself periodically over time with an angular frequency ω = 0.35 (FIG. 12a).
A variation of the elasticty number can produce interesting changes. While, a decrease in ϑ to 0.02 leads to a simple steady solution with a single transverse roll (that for the sake of brevity we omit from the discussion), vice versa, increasing it to 0.05 has just the opposite effect. A first sign of this increased complexity can be gathered from the spatio-temporal map (FIG. 11b). As a distinguishing mark with respect to earlier behavior seen in FIG. 11a, the mirror symmetry with respect to the midplane z = 5 is broken. The same information can also be inferred from the spectra reported in FIG. 12b. These relate to the velocity signals provided by probes located at different positions along the spanwise direction (points 10, 11, and 12 shown in FIG. 1, i.e. three points located in the centre of the layer and evenly spaced along z). The related frequencies are clearly different; while ω 10 = 0.28, ω 11 = ω 12 = 0.556.
A further increase in ϑ has the remarkable effect of producing smaller scale details in the spatio-temporal pattern (FIG. 11c and FIG. 11d), i.e. features that correspond to higher values of the wavenumber in space and, accordingly, higher frequencies in time (FIG. 12c  and FIG. 12d).
As a final look at FIG. 12d would indicate, moreover, signals probed at different locations follow slightly different laws, indeed, while the signal related to probe 10 Angular frequency (ω) with the higher PSD as a function of the elasticity number ϑ. The × simply marks the points in the (ϑ, ω) plane, • marks the cases that display a spectrum that is more chaotic than the majority of the other points and denotes a case where two frequencies have the same PSD.
scales as the 2D cases analyzed in the preceding section (−3.2/ − 6.4 law), the other two signals (11 and 12) do not change the inclination (they obey a ω −3.2 scaling). Most importantly, the corresponding spatiotemporal map (FIG. 12d) also illustrates that in this case the 3D disturbances have a predominant tendency to travel in the downstream direction (we will consider the important implications of this apparently cursory observation in section V). The remainder of this section is finally devoted to the analysis of the other branch of solutions, i.e. the cases with Ma = 2500. In particular, three values of ϑ are considered, namely 0.01, 0.02 and 0.05 (in this regard, we wish to recall that the companion 2D simulations carried out for ϑ = 0.01 and 0.02 revealed the trivial steady solution with a single cell occupying the entire length of the cavity, while the elastic wave traveling downstream was found for ϑ = 0.05).
Once again, as witnessed by the outcomes of the 3D simulations, consideration of the spanwise direction can cause a dramatic departure from the dynamics obtained under the constraint of 2D flow. In particular, unlike the corresponding stationary or periodic 2D counterparts, all these cases exhibit a chaotic behavior (not shown for the sake of brevity).

V. DISCUSSION
As made evident by the descriptions reported in the earlier section, the problem of transition in a viscoelastic fluid layer supporting Marangoni stresses is rich, and several issues contribute to make it more involved. In the present section an attempt is made to shed some additional light on these dynamics through an interesting analogy with 'other phenomena' and adequate consideration of the dominant theories available in the literature for the interpretation of fluid phenomena with elasticity.
In particular, in our endeavor to meet these objectives, the discussion is articulated along two different threads, namely, first, we consider the interconnections between the dimensionality of the disturbances and the observed sequence of bifurcations and, second, we analyze how their specific nature also contributes to determine the progression towards chaos.

A. Spatial and Physical Nature of the observed disturbances
In this subsection, the details of 'how' and 'when' disturbances which are intrinsically 2D or 3D in nature are excited, compete and eventually cause the onset of turbulence, are the main subject of discussion. In order to do so, we initially transcend the specific nature of the mechanisms causing instability and focus on purely spatial aspects, i.e. the different orientation of the recognizable multi-cellular structures produced in both 2D and 3D as the elasticity parameter is increased. In this regard, we rely on a 'similarity approach', that is, some arguments are developed in the light of the similar dynamics which have been observed in a companion category of phenomena, namely, the instabilities of buoyancy (gravitational convection) in differentially heated (along the horizontal direction) liquid metals. This subject, also known as the Hadley flow problem, has attracted much attention over the last decades, resulting in a variety of well-established results.
For relatively small values of the temperature difference, the Hadley flow displays a basic state that is very similar to that of Marangoni convection, i.e. a horizontally elongated single circulation extending between the hot and cold sides of the fluid container. As the control parameter is increased, it can support disturbances with transverse and/or longitudinal orientation just like those described in section IV. In other words, an effort is made here towards the identification of classes of "universality" that do not depend on the inertial or viscoelastic nature of the considered problem, which however may help to explain the present findings or provide useful hints or clues. This is the reason why in the following some relevant historical details about this specific category of buoyant flow are briefly recalled.
In this regard it is certainly worth starting from the simple remark that Hart 78 was the first to assess the response of differentially heated liquid metals to fluiddynamic disturbances having a transverse or longitudinal orientation. Later, Gill 79 examined in more detail the second category of disturbances. Taken together, these initial investigations provided relevant information about the expected patterning behavior and the mechanisms responsible for the excitation of these perturbations. It was shown that while in the first case (transverse disturbances) 2D circulations appear close to the inflection point of the basic velocity profile (such perturbation rolls therefore develop in a direction perpendicular to the basic flow), in the latter case, the axes of the rolls emerging as a result of the instability are parallel to the basic flow.
Put simply, the longitudinal rolls emerging in liquid metals combine with the basic unicellular flow typical of the Hadley convection thereby forcing the fluid parcels to describe helical paths in space. For the sake of completeness, it is also worth mentioning that, as yet shown by these landmark studies, while the instability leading to transverse rolls is driven by the mean shear stress (this is the reason why it is often referred to as "shear instability" and the related disturbances as hydrodynamic), the longitudinal rolls are made possible by the dynamical coupling between the mean shear stress and the buoyancy force, i.e. thermal effects directly contribute to the instability mechanism, from which the denomination of "hydrothermal disturbances". Subsequent instructive efforts appearing in the literature have clarified that these modes of convection are not mutually exclusive, nor are they progressive and that the development of turbulence via a hierarchy of instabilities can involve a rich variety of concurrent paths or lines of evolution. Although a plethora of studies have been conducted on these specific aspects (we apologize to all whose work is not included in this account), here we refer to the works by Lappa and Ferialdi, because in those analyses the dynamics for a small value of the Prandtl number (P r = 0.01) were investigated with and without the constraint of 2D flow, Refs. 80 and 81 respectively (thereby following an approach similar to that undertaken in the present work).
Given these premises, the sought similitude with the phenomena presented in section IV can therefore be introduced as follows: 1) in Ref. 80 2D disturbances were manifesting as waves traveling in the fluid from the hot side towards the cold one, 2) 3D disturbances spontaneously produced by the numerical simulations 81 were causing a sudden transition to relatively chaotic fully three-dimensional states for values of the controlling parameter (the Rayleigh number) smaller than those needed to excite regular (time-periodic) streamwise waves in the 2D case. The latter perturbations were clearly displaying a 3D nature as witnessed by the presence of recognizable velocity peaks along the third direc- tion and spatially extended vortices along the spanwise direction. Through numerical experiments based on the use of different initial conditions, it was concluded that in 3D the spatially pervasive presence of one mode of convection does not prevent the system from developing in parallel disturbances pertaining to the other category and vice versa, thereby leading to"hybrid states", which can deeply influence the path of progression towards chaos 81 .
Using this knowledge as an interpretation key, it is easy to recognize that the phenomena described in section IV obey very similar dynamics. value of the elasticity number close to the ϑ cr predicted through the 2D analysis is exceeded, a sudden variation can be noticed in the 3D results. This change is qualitatively and quantitatively substantiated by the increased number of recognizable features along the streamwise direction in FIG. 11c and 11d. These can be ascribed to the emergence of a train of corrugations or transverse rolls like those visible in FIG. 4c and 4d, respectively (which means that 2D disturbances have been excited and they coexist with the longitudinal ones). The increased spatial complexity of the pattern obviously stems from the helical paths that result from the superposition of transverse and longitudinal convective features. The increasing role played by the transverse (elastic) waves as the elasticity parameter grows, however, can also be inferred from the temporal behavior of the dominant disturbance, which (as explained at the end of section IV) displays an increased tendency to travel in the downstream (streamwise) direction for larger values of ϑ (FIG. 11d).
As one may infer yet by building on the companion buoyant phenomena in liquid metals, the coexistence of the two classes of disturbances in viscoelastic Marangoni flow is not limited to purely spatial effects. It has also a remarkable impact on the mechanisms that lead to turbulence; this is clearly witnessed by the smaller value of the elasticity number required in 3D to get an involved frequency spectrum (like those shown in FIG. 12c and  12d). Using arguments analogous to those elaborated by Lappa and Ferialdi 81 , a simple way to think about this trend lies in considering the non-linear interplay between perturbations that pertain to different categories.
In other words, just like the collision of two or more limit cycles can lead to the emergence of a strange attractor for the dynamics of buoyant liquid metal flow, independent disturbances (multiple solutions), having frequencies that are not commensurate, can produce a chaotic frequency spectrum in viscoelastic Marangoni flow for relatively small values of the Marangoni and elasticity numbers 19 (refer to section V B for additional arguments about the turbulence observed in the present study).
Remarkably, the classes of universality identifiable through this analogy can be pursued even further. Additional comparison of the dynamics investigated in Refs. 80 and 81 and the present ones, indeed, indicates that in both cases the 2D disturbances do not depend on the properties of the temperature field. The latter serves only as a driver to produce the physical effect required to put the fluid in motion (the buoyancy force for the liquid metals and the thermocapillary stresses in the present case). In other words, thermal effects do not play an active role in the instability mechanism.
As witnessed by the distribution of elastic energy in FIG. 13, the main physical process underpinning the instability in the present case is the interaction between the streamline curvature and the elastic stresses due to the stretching of polymer molecules.
The purely elastic nature of the hierarchy of bifurcations responsible for the onset of chaos in these cases is also indirectly demonstrated by the additional results obtained by repeating the 2D simulations after freezing the temperature field in a configuration corresponding to the basic steady state. In this regard, FIG. 14 clearly shows, that although thermal disturbances are filtered out, instabilities can still occur and produce a chaotic state with properties similar to those obtained in the fully coupled case. Remarkably, this may be regarded as an additional feature supporting the parallelism between the 2D hydrodynamic and elastic disturbances developed here.
Additional insights into the affinity between the two distinct categories of flows stem naturally from a comparison of FIG. 13 and FIG. 15.
Unlike the two-dimensional ones, the longitudinal disturbances emerging in viscoelastic Marangoni flow draw energy from thermal effects, as made evident by the spatial correlation of the elastic energy with the temperature field in the xz-plane (and with the surface stresses in the xy-plane). A more rigorous verification of this physical connection can be gained once again by decoupling the velocity and temperature fields, the reader being referred to the outcomes of the related 3D simulations shown in FIG. 16. The two-dimensional nature of the emerging flow and its structure (a single slightly corrugated cell) is instrumental in proving that, like the longitudinal modes of buoyancy convection in liquid metals require a dynamical coupling between the mean shear stress and the buoyancy force (a dynamical balance between the inertial and gravitational forces that makes the role played the thermal effects significant in the instability mechanism), the equivalent longitudinal modes in viscoelastic Marangoni flow rely on a similar coupling between the elastic and the thermocapillary forces.

B. Elastic turbulence
Still pursuing the parallelism introduced in section V A, we now concentrate on the differences, i.e on those aspects which set elastic turbulence apart from the corresponding phenomena obtained when the fluid has no elastic properties.
In particular, in order to clarify the distinguishing marks, first we consider again the case of turbulence in buoyant liquid metals (Newtonian fluids for which the concepts related to the so-called inertial turbulence apply), and then move to the viscoelastic case, still invoking the relevant literature and introducing additional levels of complexity as required.
As developed in the following, the interpretations for these two types of turbulence involve two types of relations, namely that between inertial and dissipation effects in Newtonian fluids and that between polymer molecule deformation and fluid flow in viscoelastic fluids. The former is well-known and, indeed, the underlying theory dates back to the seminal works by Kolmogorov, who theorized the existence of a hierarchy of scales through which the kinetic energy flowing in the system per unit time is balanced precisely by the amount of energy dissipated per unit time (while energy cascades at a constant dissipation rate from macroscopic phenomena towards smaller scales 82,83 ).
Using these concepts, this author elaborated a very interesting picture of turbulence by arguing that in the chaotic scale-reduction mechanism, at a certain distance from the largest scale, the macroscopic directional biases are lost allowing therefore turbulence to become homogeneous and isotropic, i.e. direction-independent. This is equivalent to stating that in the chaotic scalereduction mechanism, at a certain stage the flux of the cascading quantity across any scale becomes a function only of dynamic variables on that scale, which allows the energy spectrum to take a universal behavior (this being E(ω) ∝ ω −5/3 in the case of inertial turbulence, where the fixed exponent −5/3 accounts for the universal dependence of the energy on the "local", i.e. scaledependent, wavenumber or frequency, 84 ).
Although these concepts have been found to describe properly the known properties of turbulence on relatively small scales (in the so-called "inertial range") for several Newtonian flows (including the liquid metals discussed before and examined by Ref 81 for the case of competing transverse and longitudinal disturbances), however, they cannot be used to describe elastic turbulence in viscoelastic fluids (see, e.g., Refs. 28, 85, and 86).
As also confirmed by the present results for Marangoni flow, elastic turbulence does not depend on the Reynolds number or on equivalent non-dimensional groups. Rather chaos is enabled on increasing the elasticity (nondimensional) parameter, which in turn is directly proportional to the relaxation time λ. The latter accounts for the ability of viscoelastic stresses to survive (unlike viscous stresses) in the limit as the fluid approaches a motionless state. This is actually the mechanism by which initial disturbances can lead to the onset of turbulence. These can produce polymer molecules stretching, the deformation of the molecules can cause secondary flows which further stretch the molecules, thereby allowing the amplification of an initial small disturbance through an iterative cause-and-effect coupling mechanism 85,87 .
As implicitly revealed by these arguments, one should therefore expect the evolution of turbulence not to be determined by the assumption of local balance between small scale (in the inertial range) energy production and dissipation, rather to be somehow controlled by the alternate equilibrium between kinetic energy and elastic energy (the energy cascading from the injection scale and dissipating due to polymer relaxation). This is indeed the rationale of various theories elaborated more recently, where using these arguments it has been found (see, e.g. Ref 88 and 89) that the velocity spectrum E(ω) should decay faster than ω −3 (which is in agreement with available experiments 87 and the present results).
In order to put the present results for viscoelastic flow in perspective, it is therefore worth referring to this category of studies and, in particular, the very recent work by Gupta et al. 60 , where although the analysis was conducted under the constraint of two-dimensionality, a kind of flow that displays a remarkable affinity with Marangoni convection was considered, i.e. fluid motion inside a lid driven cavity.
For a fluid Prandtl number equal to that considered in the present work, Gupta et al. 60 observed an interesting dependence of the fitted values of the power-law exponent on the considered Reynolds number. Indeed, the exponent was found to be −3.18 in many situations. We wish also to recall that the so-called elasto-inertial turbulence is expected to display a −14/3 scaling, which is far from the Kolmogorov scaling of −5/3 typical of inertial turbulence, but relatively close to the −3.5 scaling observed for purely elastic turbulence (see, e.g., Refs. [90][91][92]. Although in the frame of the present study we have not observed turbulence satisfying the −14/3 scaling, as a concluding remark, we would like to highlight that, apparently, regions exist in the considered fluid layer where the typical behavior of inertial turbulence is recovered. In particular, as evident in FIG. 17, this happens in proximity to the free surface where the driving force is located and velocity attains its highest magnitude. As one may expect on the basis of simple ("heuristic") considerations, the interval of frequencies over which the behavior is well described by the Kolmogorov exponent grows in the situation with the larger value of the Marangoni number. In both cases, however, the Kolmogorov scaling shows up only over a distance (starting from the free surface) corresponding to approximately one quarter of the total fluid layer depth. More precisely, while for Ma = 500 the change in the exponent (from the ≈ −3.2 typical of the bulk to −5/3) occurs at 20%, for Ma = 2500 a more complex behavior is obtained with the exponent taking the values ≈ −3.2 for y ≤ 0.75, ≈ −2.5 for y = 0.8, ≈ −5/3 for y = 0.85 and ≈ −1.1 for y = 1.

VI. CONCLUSIONS
The emerging properties of Marangoni convection in a layer of viscoelastic fluid have been investigated numerically to draw some general conclusions about the behavior of this category of flows in the non-linear regime (i.e. after the disturbances have saturated their amplitude). It has been found that the elasticity of the fluid causes a remarkable decrease in the value of the Marangoni number required to excite oscillatory flow (i.e. a Hopf bifurcation) regardless of the problem dimensionality (i.e. regardless of whether the flow is constrained to retain a two-dimensional behavior or not). This is reminiscent of the corresponding "overstable" phenomena in the companion category of Marangoni-Bénard flows. The patterning behavior in the fully 3D case is driven by the interplay of two distinct categories of perturbations, which can be distinguished according to the space orientation of the emerging convective structures and the nature of the mechanisms feeding them. These can be classified accordingly as transverse or longitudinal modes of convection. The former show up as transverse waves that travel in the same direction of the imposed temperature difference (in the downstream direction, as opposed to the upstream propagation sense of classical hydrothermal waves) and owe their existence to a purely elastic instability induced by the curvature of the streamlines. By contrast, the latter rely on the coupling between the temperature field and elastic effects established through the balance of stresses at the free interface. These manifest themselves as longitudinal rolls that combine with the basic unicellular flow typical of the Marangoni convection, thereby forcing the fluid parcels to describe helical paths in space.
Building on an interesting parallelism between the typical hierarchy of bifurcations of buoyant convection (the Hadley flow) in liquid metals and that of Marangoni flow in viscoelastic liquids, it has been shown that in both cases, on increasing the magnitude of the driving force, longitudinal disturbances are selected first and then transverse perturbations are excited and coexist with the longitudinal ones. This results in a set of independent convective modes that can produced turbulence through non-linear interaction.
Entering the turbulent regime, however, does not require necessarily an increase in the intensity of the driving force (the Marangoni number). For a fixed magnitude of the driving force, chaos can also be excited by increasing the level of elasticity. Accordingly, the emerging non-linear states display the typical properties of elas-tic turbulence as witnessed by the scaling laws and the related exponent ≈ −3.2 recognizable in the frequency spectrum (as opposed to the −5/3 exponent typical of inertial turbulence).
As shown by the present simulations, however, these two forms of turbulence are not mutually exclusive. These can coexists and cause a gradual transition in the aforementioned exponent along a direction perpendicular to the free interface, which indicates that the dichotomy often drawn between the related physical mechanisms should not be taken in a strict sense when dealing with viscoelastic Marangoni flow in non creeping conditions. An exciting prospect for the future is the generalization of these findings to other configurations, which have enjoyed a widespread use in the literature for the investigation of the fundamental properties of Marangoni convection, such the liquid bridge or the classical annular pool used to mimic the CZ crystal growth process.