2 Formulation

We consider the evolution of an incompressible fluid inside a vertical right circular cylindrical container of diameter $D$ and height $H$. The lower horizontal plate is maintained at a higher temperature than the upper plate, while the circular walls are assumed to be thermally insulating. All walls are impenetrable and no-slip. The layer rotates anticlockwise about the $z$ axis, pointing vertically upwards, with a constant angular velocity $\unicode[STIX]{x1D734}=\unicode[STIX]{x1D6FA}\boldsymbol{e}_{z}$, while gravity points downwards, $\boldsymbol{g}=-g\boldsymbol{e}_{z}$. The kinematic viscosity $\unicode[STIX]{x1D708}$ and thermal diffusivity $\unicode[STIX]{x1D705}$ are assumed to be constant.

In the Boussinesq approximation, using the rotation time $1/(2\unicode[STIX]{x1D6FA})$ as a unit of time and the depth $H$ of the layer as a unit of length, the dimensionless equations are

(2.1)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\unicode[STIX]{x1D735}p-\boldsymbol{e}_{z}\times \boldsymbol{u}+\frac{RaE^{2}}{Pr}\unicode[STIX]{x1D703}\boldsymbol{e}_{z}+E\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

(2.2)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(2.3)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}=u_{z}+\frac{E}{Pr}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}. & \displaystyle\end{eqnarray}$$

Here $\boldsymbol{u}=(u_{x},u_{y},u_{z})$ is the velocity, $p$ is the pressure and $\unicode[STIX]{x1D703}$ is the temperature fluctuation relative to the conduction profile, which varies linearly with respect to the vertical coordinate $z$. The control parameters are the Rayleigh number $Ra=\unicode[STIX]{x1D6FC}g\unicode[STIX]{x0394}TH^{3}/(\unicode[STIX]{x1D708}\unicode[STIX]{x1D705})$, the Ekman number $E=\unicode[STIX]{x1D708}/(2\unicode[STIX]{x1D6FA}H^{2})$ and the Prandtl number $Pr=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}$. These dimensionless quantities involve $\unicode[STIX]{x1D6FC}$, the coefficient of thermal expansion, and $\unicode[STIX]{x0394}T$, the imposed temperature difference between the two horizontal plates. The aspect ratio $\unicode[STIX]{x1D6E4}\equiv D/H$ specifies the geometry.

Figure 1.

Summary of the runs performed for this study (full symbols, red circle), along with past numerical (Scheel et al. Reference Scheel, Paul, Cross and Fischer2003; Sánchez-Álvarez et al. Reference Sánchez-Álvarez, Serre, Crespo del Arco and Busse2005; Horn & Schmid Reference Horn and Schmid2017; de Wit et al. Reference de Wit, Aguirre Guzman, Madonia, Cheng, Clercx and Kunnen2020; Zhang et al. Reference Zhang, van Gils, Horn, Wedi, Zwirner, Ahlers, Ecke, Weiss, Bodenschatz and Shishkina2020) and experimental (Ecke et al. Reference Ecke, Zhong and Knobloch1992; Aurnou et al. Reference Aurnou, Bertin, Grannan, Horn and Vogt2018; de Wit et al. Reference de Wit, Aguirre Guzman, Madonia, Cheng, Clercx and Kunnen2020) studies. The lower black line corresponds to the onset of wall modes as predicted by linear theory, while the upper one corresponds to the onset of bulk modes. Our study focuses mostly on the grey region in between, where only the wall mode is unstable. Note that the studies reported on this plot (symbols) do not all have the same Prandtl number $Pr$ nor the same aspect ratio $\unicode[STIX]{x1D6E4}$.

We solve equations (2.1)–(2.3) using the open-source spectral-element code Nek5000 (Fischer Reference Fischer1997) (NEK5000 Version 19.0, Argonne National Laboratory, Illinois; available at https://nek5000.mcs.anl.gov). The mesh is composed of up to 19 968 hexahedral elements and we use a polynomial order up to $N=16$ including dealiasing. The mesh is refined close to the boundaries of the domain in order to properly resolve both Ekman and Stewartson layers. Numerical convergence of the results has been checked by increasing the polynomial order.

3 Nonlinear dynamics of wall modes

We are interested in the dynamics of the wall modes, which are the first unstable modes in laterally bounded rotating RBC. The critical Rayleigh number for the appearance of wall modes can be found in Herrmann & Busse (Reference Herrmann and Busse1993) and Kuo & Cross (Reference Kuo and Cross1993) (see Zhang & Liao (Reference Zhang and Liao2009) for higher-order corrections) and is given in our dimensionless units by $Ra_{c}^{wall}\approx \unicode[STIX]{x03C0}^{2}(6\sqrt{3})^{1/2}E^{-1}$, while the critical Rayleigh number for the onset of the bulk mode is given by $Ra_{c}^{bulk}\approx 3(\unicode[STIX]{x03C0}^{2}/2)^{2/3}E^{-4/3}$ (Chandrasekhar Reference Chandrasekhar1961; Clune & Knobloch Reference Clune and Knobloch1993; Liao, Zhang & Chang Reference Liao, Zhang and Chang2006). These expressions are valid for low Ekman numbers and are independent of both the Prandtl number and aspect ratio provided the latter is sufficiently large. Although formally obtained for different boundary conditions at the top and bottom, the analysis of Clune & Knobloch (Reference Clune and Knobloch1993) explains why the leading-order behaviour is independent of these boundary conditions. Both expressions are shown in figure 1, where we also report previous experimental and numerical studies of interest. Interestingly, the range of Rayleigh numbers for which only wall modes are unstable increases as the Ekman number decreases. The above asymptotic results belie the fact that the onset azimuthal wavenumber is $m=(\unicode[STIX]{x03C0}/2)\unicode[STIX]{x1D6E4}(2+\sqrt{3})^{1/2}\approx 3\unicode[STIX]{x1D6E4}$ as $E\rightarrow 0$ (Zhang & Liao Reference Zhang and Liao2009). Since $m$ is necessarily an integer, each $m$ defines a neutral stability curve, and the intersection between two successive curves heralds a transition from onset wavenumber $m$ to onset wavenumber $m+1$ as $\unicode[STIX]{x1D6E4}$ increases. Thus the critical Rayleigh number is in fact an oscillatory function of $\unicode[STIX]{x1D6E4}$, although its minima are independent of $\unicode[STIX]{x1D6E4}$ for large $\unicode[STIX]{x1D6E4}$.

Figure 2.

Left: Visualizations of the instantaneous vertical velocity $u_{z}$ (left side) and fluctuating temperature $\unicode[STIX]{x1D703}$ (right side) in the mid-plane $z=0.5$ in a $\unicode[STIX]{x1D6E4}=1.5$ cylinder in the quasi-stationary state. Middle: Instantaneous temperature fluctuations $\unicode[STIX]{x1D703}$ on the sidewall $r=\unicode[STIX]{x1D6E4}/2$ as a function of the azimuthal angle $\unicode[STIX]{x1D719}$ and vertical coordinate $z$ in the quasi-stationary state. Right: Spatio-temporal plots showing the evolution of the temperature fluctuation $\unicode[STIX]{x1D703}$ at $z=0.5$ and radius $r=0.74$, with $\unicode[STIX]{x1D719}$ plotted horizontally and time $t$ (in rotation units) increasing downwards. The time interval is different in each panel. In each case, the simulation is run for at least 100 free-fall times, which corresponds to thousands of rotation times. The parameters are $\unicode[STIX]{x1D6E4}=1.5$, $E=10^{-6}$ and $Pr=1$; in such a cylinder, the wall is at $r=0.75$. The Rayleigh number increases from top to bottom: (a) $Ra=5\times 10^{7}$, (b) $Ra=2\times 10^{8}$, (c) $Ra=5\times 10^{8}$ and (d) $Ra=2\times 10^{9}$. In (b), the dashed horizontal line indicates a reduction of $Ra$ from $2\times 10^{8}$ to $5\times 10^{7}$: the $m=4$ mode survives, albeit with a slower precession frequency.

In this section, we fix $E=10^{-6}$ so that $Ra_{c}^{wall}\approx 3.2\times 10^{7}$ while $Ra_{c}^{bulk}\approx 8.7\times 10^{8}$. We then vary the Rayleigh number from $Ra=5\times 10^{7}$ to $Ra=8\times 10^{9}$ to explore the onset of the wall mode, its nonlinear saturation and finally its coexistence with the bulk mode. The aspect ratio is $\unicode[STIX]{x1D6E4}=1.5$ and the Prandtl number is $Pr=1$ for simplicity (see § 4 for parameters relevant to the recent experiments). The simulations are initialized with small-amplitude temperature fluctuations. In each case we observe exponential growth followed by nonlinear saturation. Visualizations of the vertical velocity and the temperature fluctuation in the mid-plane $z=0.5$ in the saturated state are shown in the left panels of figure 2. At $Ra=5\times 10^{7}$, i.e. relatively close to the onset of the wall mode, we find a retrograde $m=6$ travelling wave localized near the sidewall (figure 2a). The travelling nature of the wall mode is reflected in the asymmetry of its instantaneous profile (see the temperature fluctuations on the side boundary in the middle panels of figure 2), but is best seen by looking at spatio-temporal diagrams of the temperature fluctuation in the mid-plane $z=0.5$ and an arbitrary location $r$ close to the boundary: $r=0.74<\unicode[STIX]{x1D6E4}/2=0.75$ (right panels in figure 2). The results at other radii and depths are very similar provided one remains sufficiently close to the outer boundary. This travelling mode is similar to that observed in experiments in a $\unicode[STIX]{x1D6E4}=2$ cylinder rotating at lower rotation rates (Zhong et al. Reference Zhong, Ecke and Steinberg1991, Reference Zhong, Ecke and Steinberg1993). The structure of this mode is also similar to that of the critical mode for $\unicode[STIX]{x1D6E4}=4$, $E=5\times 10^{-5}$ and $Pr=1$ computed by Zhang & Liao (Reference Zhang and Liao2009) for which $m=10$ while the asymptotic wavenumber as $E\rightarrow 0$ is $m\approx 12$. In contrast, in our case, computed for $E=10^{-6}$, the observed wavenumber is larger than the onset theoretical wavenumber $m\approx 4.5$. Evidently our choice of parameters is outside the asymptotic regime $E\rightarrow 0$ and $\unicode[STIX]{x1D6E4}\rightarrow \infty$.

As the Rayleigh number increases – while remaining below the threshold for the bulk mode – the wall mode becomes fully nonlinear and more complex dynamical behaviour is seen. We first observe (see figure 2b) an Eckhaus-type instability (Janiaud et al. Reference Janiaud, Pumir, Bensimon, Croquette, Richter and Kramer1992; Knobloch Reference Knobloch1994), as found experimentally at much higher Ekman number by Liu & Ecke (Reference Liu and Ecke1999). This instability leads to the gradual coarsening or merging of adjacent cells but the flow remains laminar. The transitions between different azimuthal wavenumbers are expected to be strongly hysteretic for reasons explained in Knobloch (Reference Knobloch1994), a behaviour that is illustrated in figure 2(b), where the $m=4$ wall mode persists even when the Rayleigh number is reduced from $2\times 10^{8}$ back to $5\times 10^{7}$, i.e. to a Rayleigh number at which $m=6$ is the first mode that appears. At even higher Rayleigh numbers ($Ra=5\times 10^{8}$, see figure 2c), the $m=4$ travelling wave pattern persists but becomes fully nonlinear and starts to emit plumes into the stable bulk resulting in short-wavelength perturbations of the mode. These fluctuations tend to propagate in the opposite direction from that of the travelling wave, i.e. in a cyclonic sense. Finally, a $m=3$ wall mode survives even for Rayleigh numbers beyond the onset of bulk convection, where it takes the form of a nonlinear wall-confined travelling wave superposed on small-scale bulk turbulence (figure 2d).

From the spatio-temporal diagrams one can easily extract the drift frequency $\unicode[STIX]{x1D714}_{d}$ as a function of the Rayleigh number. To do so, we perform a temporal Fourier transform of the fluctuating temperature at a given point and extract the frequency corresponding to the maximum of the power spectrum. This frequency is shown in figure 3(a) (all values are negative so we plot $-\unicode[STIX]{x1D714}_{d}$). The drift frequency increases continuously with increasing $Ra$ and connects to the theoretical prediction of Herrmann & Busse (Reference Herrmann and Busse1993) at onset, $\unicode[STIX]{x1D714}_{c}\approx -59E/Pr$. This increase appears to follow a linear scaling with the Rayleigh number, $|\unicode[STIX]{x1D714}_{d}-\unicode[STIX]{x1D714}_{c}|\sim Ra-Ra_{c}^{wall}$, as observed at higher Ekman numbers by Horn & Schmid (Reference Horn and Schmid2017) and more recently by de Wit et al. (Reference de Wit, Aguirre Guzman, Madonia, Cheng, Clercx and Kunnen2020) (see § 4 below). In fact such a linear relationship between the wall mode drift frequency and the Rayleigh number is expected close to the Hopf bifurcation from which the mode originates (Ecke et al. Reference Ecke, Zhong and Knobloch1992; Goldstein et al. Reference Goldstein, Knobloch, Mercader and Net1993), but it is unexpected that it continues to hold over nearly two decades in Rayleigh number, and extends into the regime where the bulk mode is unstable and even turbulent. This is a first indication of the robustness of the wall state.

All wall states can be decomposed into an axisymmetric component with wavenumber $m=0$ and a drifting non-axisymmetric component with wavenumbers $m>0$. We identify the former with the observed zonal flow (Sánchez-Álvarez et al. Reference Sánchez-Álvarez, Serre, Crespo del Arco and Busse2005; de Wit et al. Reference de Wit, Aguirre Guzman, Madonia, Cheng, Clercx and Kunnen2020; Zhang et al. Reference Zhang, van Gils, Horn, Wedi, Zwirner, Ahlers, Ecke, Weiss, Bodenschatz and Shishkina2020), and use the latter to identify the drift frequency $\unicode[STIX]{x1D714}_{d}$. To study the zonal flow, we compute the time and volume average of the azimuthal velocity component $u_{\unicode[STIX]{x1D719}}$ in the quasi-stationary state (restricting the volume average to a region near the sidewall leads to qualitatively similar results). As shown in figure 3(b), the mean zonal flow amplitude also increases linearly with $Ra-Ra_{c}^{wall}$ and becomes vanishingly small close to onset. This observation is consistent with the fact that the Nusselt number also increases linearly with the supercriticality parameter near onset (not shown; see Zhong et al. (Reference Zhong, Ecke and Steinberg1993) and Horn & Schmid (Reference Horn and Schmid2017)). Note that the linear scaling once again extends over the whole range of Rayleigh numbers considered here, including cases with and without bulk turbulence, with no significant impact of the transition between them. This scaling is, of course, not expected to hold indefinitely as the Rayleigh number is increased even further; the flow becomes less rotationally constrained and contributions from the bulk eventually dominate.

Figure 3.

(a) Drift frequency $-\unicode[STIX]{x1D714}_{d}$ as a function of $Ra$ for $\unicode[STIX]{x1D6E4}=1.5$, $E=10^{-6}$ and $Pr=1$. The results for the full cylinder (▫) and the cylinder with a barrier (▵, see § 5 below) coincide. The theoretical value $\unicode[STIX]{x1D714}_{c}\approx -59E/Pr$ predicted by Herrmann & Busse (Reference Herrmann and Busse1993) for the onset of the instability in the presence of a planar wall is also reported (open circle). The dot-dashed line corresponds to the linear scaling $|\unicode[STIX]{x1D714}_{d}|=|\unicode[STIX]{x1D714}_{c}|+c(Ra-Ra_{c}^{wall})$ with $c$ an arbitrary constant. (b) Corresponding volume- and time-averaged zonal velocity $\langle u_{\unicode[STIX]{x1D719}}\rangle _{V,t}$. The dot-dashed line corresponds to the scaling $\langle u_{\unicode[STIX]{x1D719}}\rangle _{V,t}\sim Ra-Ra_{c}^{wall}$. The volume integration is performed only over the left half of the cylinder for the cases with barrier (see § 5). (c) Azimuthally, vertically and temporally averaged zonal velocity $\langle u_{\unicode[STIX]{x1D719}}\rangle _{\unicode[STIX]{x1D719},z,t}$ as a function of the radial coordinate $r$. Positive values correspond to cyclonic motions while negative values correspond to anticyclonic motions. The two vertical lines indicate the Stewartson layer scales $E^{1/3}$ and $E^{1/4}$ (Stewartson Reference Stewartson1957). (d) Divergence of the Reynolds stress in the azimuthal direction for the case $Ra=5\times 10^{8}$. Results are computed at a given time during the exponential phase of the instability (arbitrarily rescaled) and time-averaged during the nonlinear saturated phase.

Following de Wit et al. (Reference de Wit, Aguirre Guzman, Madonia, Cheng, Clercx and Kunnen2020), we determine the radial structure of the zonal flow by computing the azimuthally, vertically and temporally averaged zonal velocity $\langle u_{\unicode[STIX]{x1D719}}\rangle _{\unicode[STIX]{x1D719},z,t}$ as shown in figure 3(c). Close to onset, the mean azimuthal flow is negligible, but with increasing $Ra$ it strengthens rapidly and develops a concentric ring structure, with cyclonic flow in the outer ring and anticyclonic flow in the inner ring. This flow structure is typical of Stewartson layers (Stewartson Reference Stewartson1957; Kunnen et al. Reference Kunnen, Stevens, Overkamp, Sun, van Heijst and Clercx2011), with the cyclonic mean flow confined to a layer of thickness $E^{1/3}$ close to the boundary and the anticyclonic flow present in a broader layer of thickness $E^{1/4}$, exactly as observed by Zhang et al. (Reference Zhang, van Gils, Horn, Wedi, Zwirner, Ahlers, Ecke, Weiss, Bodenschatz and Shishkina2020) (see their figure 3b). Interestingly, one can now compute the typical Reynolds number associated with this mean differential rotation. Based on the largest layer thickness $E^{1/4}$ and the azimuthal mean velocity difference in the profile shown in figure 3(b), one obtains $Re\approx 2.3$ at $Ra=2\times 10^{8}$ and $Re\approx 9.7$ at $Ra=5\times 10^{8}$. Since the critical Reynolds number for shear barotropic instabilities is typically around $Re\approx 10{-}20$ independently of $E$ (Hide & Titman Reference Hide and Titman1967; Niino & Misawa Reference Niino and Misawa1984; Früh & Read Reference Früh and Read1999), we conjecture that the small-scale fluctuations observed at $Ra=5\times 10^{8}$ but not at $Ra=2\times 10^{8}$ are driven by a shear instability of the mean barotropic zonal flow (Busse Reference Busse1968), although a complete analysis of this instability is beyond the scope of this paper. The same nested cyclonic–anticyclonic structure is observed even after the appearance of the bulk state, indicating that this structure is a consequence of the survival of the nonlinear wall state into this regime. The observed spreading of the anticyclonic flow into the bulk as $Ra$ increases, well beyond the linear Stewartson layer, indicates that additional effects come into play. This spreading is observed as soon as the wall mode becomes unstable to secondary instabilities (i.e. for $Ra\gtrsim 5\times 10^{8}$), even before the bulk mode sets in, and is therefore a property of the nonlinear wall state that cannot be fully attributed to the emergence of the bulk state.

The mean zonal flow is driven by the non-axisymmetric structure of the wall state itself. While the drift frequency $\unicode[STIX]{x1D714}_{d}$ is a consequence of a fundamentally linear mechanism, the zonal flow is driven nonlinearly by Reynolds stresses. These stresses are known to sustain both zonal geostrophic and meridional non-geostrophic flows within a weakly nonlinear analysis (Liao, Zhang & Chang Reference Liao, Zhang and Chang2005). We confirm here that this mechanism continues to operate at larger Rayleigh numbers by computing the divergence of the Reynolds stresses in the azimuthally and vertically averaged zonal equation, ${\mathcal{R}}(r)=-\boldsymbol{e}_{\unicode[STIX]{x1D719}}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\langle \boldsymbol{u}^{\prime }\boldsymbol{u}^{\prime }\rangle _{\unicode[STIX]{x1D719},z})$, where $\boldsymbol{u}^{\prime }=\boldsymbol{u}-\langle \boldsymbol{u}\rangle _{\unicode[STIX]{x1D719},z}$. The results are shown in figure 3(d) for $Ra=5\times 10^{8}$. We see that the Reynolds stress associated with the linear wall mode (i.e. calculated during the exponential phase of the instability) is indeed responsible for the cyclonic mean flow within the $E^{1/3}$ Stewartson layer and for the dominant outer anticyclonic zonal flow. The structure of this stress persists into the saturated phase, even though the wall mode is now fully nonlinear. In this quasi-steady state the Reynolds stress is in balance with viscous dissipation. It remains to be seen how the stresses generated by the wall state interact with the bulk state, potentially explaining the expansion of the anticyclonic mean flow observed in figure 3(c) and in Zhang et al. (Reference Zhang, van Gils, Horn, Wedi, Zwirner, Ahlers, Ecke, Weiss, Bodenschatz and Shishkina2020).

We discuss, finally, the full spatial structure of the mean flow, which is not fully captured by the vertically averaged quantity $\langle u_{\unicode[STIX]{x1D719}}\rangle _{\unicode[STIX]{x1D719},z}$. The azimuthal and temporal average of all three components of the velocity are shown in figure 4. The intricate mean flow structure observed for $Ra=2\times 10^{8}$ persists even after the instability of the bulk mode, at $Ra=2\times 10^{9}$. Although the signal is now perturbed by the small-scale fluctuations driven by the bulk state, the same flow structure is easily recognizable, showing that the mean flow is predominantly driven by the wall mode itself even when bulk turbulence is present. As in figure 3(c), the mean flow pattern spreads out radially with increasing $Ra$, in agreement with Zhang et al. (Reference Zhang, van Gils, Horn, Wedi, Zwirner, Ahlers, Ecke, Weiss, Bodenschatz and Shishkina2020), and all components of the flow follow the linear scaling with the distance from threshold identified in figure 3(b). Moreover, complex meridional structures are present even when the zonal component dominates the mean flow. For example, the cyclonic motion observed in figure 3(c) is centred near the mid-plane and also spreads inwards as $Ra$ increases. We also recover the meridional circulation predicted by weakly nonlinear theory (Liao et al. Reference Liao, Zhang and Chang2005) with two counter-rotating cells on either side of the mid-plane. Finally, we observe strong outward radial flow inside the top and bottom Ekman layers, a direct consequence of Ekman pumping due to the differential rotation above and below.

4 Link with recent experiments and simulations at high $Ra$

Figure 4.

Temporal and azimuthal averages of all three velocity components in cylindrical coordinates as a function of the radial and vertical positions, for (a) $Ra=2\times 10^{8}$ and (b) $Ra=2\times 10^{9}$. Plots focus on the region close to the sidewall (located at $r=\unicode[STIX]{x1D6E4}/2=0.75$) and the vertical coordinate is stretched to zoom in to the upper and lower walls (the mapping $z\rightarrow \sin ^{2}(\unicode[STIX]{x03C0}z/2)$ is applied twice). The horizontal dotted lines show the vertical thickness $E^{1/2}$ of the Ekman boundary layers, while the vertical dashed lines correspond to the Stewartson layer scaling $E^{1/4}$ and $E^{1/3}$. Colour scales are chosen consistently with the linear scaling observed in figure 3(b).

Recent experiments and simulations have shown the emergence of a boundary zonal flow in confined rotating RBC (Kunnen et al. Reference Kunnen, Stevens, Overkamp, Sun, van Heijst and Clercx2011; de Wit et al. Reference de Wit, Aguirre Guzman, Madonia, Cheng, Clercx and Kunnen2020; Zhang et al. Reference Zhang, van Gils, Horn, Wedi, Zwirner, Ahlers, Ecke, Weiss, Bodenschatz and Shishkina2020). However, a link between this flow and the wall modes was not explicitly made. In further support of this link, we repeated the simulation of de Wit et al. (Reference de Wit, Aguirre Guzman, Madonia, Cheng, Clercx and Kunnen2020) using the parameters $E=10^{-7}$, $Pr=5.2$ and $\unicode[STIX]{x1D6E4}=0.2$. Instead of focusing on the dynamics of the bulk mode, as in de Wit et al. (Reference de Wit, Aguirre Guzman, Madonia, Cheng, Clercx and Kunnen2020), we consider Rayleigh numbers below the onset of the bulk mode (i.e. for $Ra<Ra_{c}^{bulk}\approx 2\times 10^{10}$), but above the onset of the wall mode ($Ra>Ra_{c}^{wall}\approx 3\times 10^{8}$). We show in figure 5(a) the vertical velocity in the mid-plane at $Ra=2\times 10^{9}$, well beyond the onset of wall modes, once the system has reached a statistically steady state. We observe an $m=1$ travelling wave, presumably a consequence of the small value of $\unicode[STIX]{x1D6E4}$ (the asymptotic onset wavenumber for $\unicode[STIX]{x1D6E4}=0.2$ is $m\approx 0.6$). This state is very similar to the boundary zonal flow observed by de Wit et al. (Reference de Wit, Aguirre Guzman, Madonia, Cheng, Clercx and Kunnen2020) (see their figure 2b, for example) at the much higher Rayleigh number of $Ra=5\times 10^{10}$. Here the nonlinear dynamics of the wall state are clearly constrained by the small aspect ratio used in these simulations compared to the cases discussed in the previous section. However, once the Rayleigh number is increased to $Ra=5\times 10^{10}$, we recover the emergence of the bulk mode superposed on the nonlinear wall mode (see figure 5b, which can be directly compared with figure 2b of de Wit et al. (Reference de Wit, Aguirre Guzman, Madonia, Cheng, Clercx and Kunnen2020)). Note that de Wit et al. (Reference de Wit, Aguirre Guzman, Madonia, Cheng, Clercx and Kunnen2020) also measured the drift frequency as a function of $Ra$ and found the following fit: $\unicode[STIX]{x1D714}_{d}\approx -6\times 10^{-10}Ra^{1.16\pm 0.06}E$ (using our rotation units).

Figure 5.

Vertical velocity in the mid-plane $z=0.5$ in a $\unicode[STIX]{x1D6E4}=0.2$ cylinder rotating at $E=10^{-7}$ with $Pr=5.2$ and (a) $Ra=2\times 10^{9}$ and (b) $Ra=5\times 10^{10}$ (this latter case can be compared with figure 2b of de Wit et al. (Reference de Wit, Aguirre Guzman, Madonia, Cheng, Clercx and Kunnen2020) at exactly the same parameters). (c) Drift frequency $-\unicode[STIX]{x1D714}_{d}$ as a function of $Ra$. The theoretical value $\unicode[STIX]{x1D714}_{d}=\unicode[STIX]{x1D714}_{c}\approx -63E/Pr$ predicted by Herrmann & Busse (Reference Herrmann and Busse1993) for the onset of the instability in the presence of a planar wall is also reported (open circle). The oblique line corresponds to the scaling $\unicode[STIX]{x1D714}_{d}\approx -6\times 10^{-10}Ra^{1.16}E$ reported in de Wit et al. (Reference de Wit, Aguirre Guzman, Madonia, Cheng, Clercx and Kunnen2020) and obtained at much larger $Ra$.

Using the same approach as in § 3, we measured the drift frequency in our simulations below the onset of the bulk mode, and in figure 5(c) we compare the results with the fit of de Wit et al. (Reference de Wit, Aguirre Guzman, Madonia, Cheng, Clercx and Kunnen2020). The results clearly bridge the gap between the onset theoretical value derived by Herrmann & Busse (Reference Herrmann and Busse1993) and the highly turbulent scaling found by de Wit et al. (Reference de Wit, Aguirre Guzman, Madonia, Cheng, Clercx and Kunnen2020), thereby providing additional evidence for the robustness of the linear scaling discussed in § 3. Note also that we reproduce many of the results presented in these studies (such as the bimodal temperature distribution and the boundary zonal flow) by considering cases below the onset of the bulk mode. This provides a clear indication that the boundary dynamics observed in these studies are related to the nonlinear saturation of the wall mode and so are not directly driven by the dynamics in the bulk (which is stable in two of our cases). The emergence of the boundary zonal flow is thus best understood as the nonlinear development of linearly unstable wall modes that eventually couple to the bulk dynamics in the manner described in these studies (Kunnen et al. Reference Kunnen, Stevens, Overkamp, Sun, van Heijst and Clercx2011; de Wit et al. Reference de Wit, Aguirre Guzman, Madonia, Cheng, Clercx and Kunnen2020; Zhang et al. Reference Zhang, van Gils, Horn, Wedi, Zwirner, Ahlers, Ecke, Weiss, Bodenschatz and Shishkina2020).

5 Robustness of the wall modes to geometry changes

The time-dependent wall modes we have described thus far are linked to the presence of a sidewall. These states behave like topologically protected states and, as shown above, survive even when the Rayleigh number is very far from its critical value, i.e. in the strongly nonlinear regime. The robustness of these states can also be demonstrated by perturbing the geometry of the domain. This has already been done in other similar systems (Dasbiswas, Mandadapu & Vaikuntanathan Reference Dasbiswas, Mandadapu and Vaikuntanathan2018; Souslov et al. Reference Souslov, Dasbiswas, Fruchart, Vaikuntanathan and Vitelli2019) and was in fact suggested in the very early studies of wall modes in rotating convection by Liu & Ecke (Reference Liu and Ecke1999); see Vasil, Brummell & Julien (Reference Vasil, Brummell and Julien2008) for a related study of drifting states in a domain of square cross-section.

Figure 6.

Vertical velocity in the mid-plane $z=0.5$ for a cylinder with a barrier. The Rayleigh number increases from (a) to (c): (a) $Ra=5\times 10^{7}$, (b) $Ra=5\times 10^{8}$ and (c) $Ra=2\times 10^{9}$. Parameters are $\unicode[STIX]{x1D6E4}=1.5$, $E=10^{-6}$ and $Pr=1$.

In this section, we return to the parameters used in § 3. However, the cylinder is now truncated and we introduce a rectangular barrier of length $\unicode[STIX]{x1D6E4}/2$ and thickness $\unicode[STIX]{x1D6E4}/20$ along a diameter with one end attached to the circular sidewall (see figure 6). Like the rest of the sidewall, the barrier is no-slip and insulating. The barrier breaks the rotation invariance of the system and introduces four discontinuities in the curvature of the sidewall. One might therefore expect that this change in the geometry of the domain will result in a drastic change in the nature of the wall state. This is not the case. As shown in figures 6 and 7, for $Ra=5\times 10^{7}$, $Ra=5\times 10^{8}$ and $Ra=2\times 10^{9}$, the wall state persists, in both its weakly nonlinear and fully nonlinear forms. More surprisingly, the drift frequency is the same as in the case without the barrier for all three Rayleigh numbers, as shown in figure 3(a). This is even more surprising when one notices that the wavelength of the wall state is comparable with the length of the barrier and larger than its thickness. As seen in figure 6, the wall mode simply wraps itself around the obstacle, following the same retrograde drift as in the case without a barrier, although its wavelength and frequency along the barrier are both slightly larger than around the rest of the cylinder (see figure 7a). Moreover, the presence of the corners appears to be responsible for increased shedding of cyclonic disturbances. Despite this, the zonal flow discussed earlier is barely modified by the presence of the barrier (see figure 7b), and its average in a direction tangential to the boundary exhibits the same nested cyclonic–anticyclonic structure along all boundaries. Nevertheless, the amplitude of the zonal flow is very slightly reduced, a fact confirmed in figure 3(b) by a reduction in its volume average, here performed over the left half of the cylinder only. Finally, as in the case without a barrier, the wall mode survives the onset of the bulk mode and remains visible in instantaneous visualizations of the vertical velocity at the mid-plane even in the presence of bulk turbulence (figure 6c).

Figure 7.

Spatio-temporal plots showing the temperature fluctuation $\unicode[STIX]{x1D703}$ at $z=0.5$ and a fixed distance $\unicode[STIX]{x1D6FF}=10^{-2}$ from the boundary at (a) $Ra=5\times 10^{7}$ and (b) $Ra=5\times 10^{8}$. The horizontal axis represents the arclength $s$ along the boundary while the vertical axis is time $t$ (in rotation units). The dotted lines indicate the positions of the four corners of the barrier. Parameters are $\unicode[STIX]{x1D6E4}=1.5$, $E=10^{-6}$ and $Pr=1$. (c) Vertical and temporal average of the velocity component tangential to the boundary $u_{\unicode[STIX]{x1D70F}}$ at $Ra=5\times 10^{8}$. We distinguish between the cylindrical boundary and the different faces of the barrier. The results are compared to the case without barrier (see figure 3c).