Flow transitions resembling bifurcations of the logistic map in simulations of the baroclinic rotating annulus

We present evidence for a sequence of bifurcations in simulations of the differentially heated baroclinic rotating annulus, similar to bifurcations associated with the logistic map. The Met. Office / Oxford Rotating Annulus Laboratory Simulation (MORALS) code is used to construct a detailed numerical regime diagram for the annulus, and the distribution of regimes in parameter space is discussed. The bifurcations are observed in a sequence of runs at high temperature forcing, identified by Poincaré sections of the dominant temperature mode amplitude time series. Higher order return maps and predictions using quadratic fits to the data are used to verify this result, and Lyapunov exponents are calculated to identify and quantify the chaotic parts of the sequence.


Introduction
The differentially heated rotating annulus [1][2][3] is a laboratory experiment that has been used for over 50 years to produce behaviour qualitatively similar to the midlatitudes of a generic planetary atmosphere. The experiment consists of two concentric cylinders maintained at different temperatures, rotating about a vertical axis, with fluid between them (Fig. 1). This setup is firmly established as an insightful laboratory analogue for certain kinds of atmospheric dynamical behaviour, and is also a useful testbed for the methods used to study them.

ACCEPTED MANUSCRIPT
Baroclinic instability, which is one of the primary mechanisms governing large scale fluid motion in the atmosphere, can be reproduced in the annulus under certain conditions when a temperature gradient is combined with rotation [1], and the laboratory setting allows this to be studied in a controlled and reproducible environment. As well as similarities with the atmosphere, the annulus displays a rich and diverse range of behaviour worthy of study in its own right, including low-dimensional dynamics. Although using the annulus as a direct analogue of the highly turbulent Earth's atmosphere may be questionable, the complex nature of real atmospheres means we need simpler analogues for testing ideas and methods. Weakly turbulent behaviour observed in the annulus is more relevant for simpler planetary atmospheres such as on Mars, where the large-scale flow is believed to be more predictable than large-scale flow on Earth [4].
This work is exclusively a computational study. Initially, we summarize the construction of a detailed numerical regime diagram for the rotating annulus; this is directly comparable with equivalent diagrams produced from earlier laboratory experiments. We then report the discovery of a previously unseen regime, whose behaviour is similar to the bifurcations produced by iteration of the logistic map [5][6][7]. In the laboratory, this type of behaviour is well documented in Rayleigh-Bénard convection cells using water and mercury: these experiments found period-doubling bifurcations in velocity and temperature time series as a result of varying external forcing, either by a temperature gradient (water) [8,9] or magnetic forcing (mercury) [10,11]. Period-doubling of a periodic amplitude modulation has also been seen as a result of baroclinic instability in a two layer, open cylinder experiment driven by differential rotation [12]. It has not been observed to date in the laboratory annulus under A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT radial thermal forcing, however. Section 2 describes the simulation. In Section 3 the regime diagram and its main features are outlined. We present evidence for a sequence of bifurcations similar to the logistic map in Section 4, and supplementary evidence using Lyapunov exponents in Section 5. Section 6 provides discussion and conclusions.

The simulation
We use the Met. Office / Oxford Rotating Annulus Laboratory Simulation (MORALS) code [13][14][15]. This solves the full Navier-Stokes, continuity, and heat transfer equations along with equations of state for density, viscosity and thermal diffusivity (Appendix A), under the Boussinesq approximation for a cylindrical annulus rotating at angular velocity Ω. The equations are cast in velocity-temperature-pressure form: u (meridional), v (zonal) and w (vertical) velocities / cm s −1 , temperature T / • C and a normalized pressure Π = P/ρ 0 / cm 2 s −2 . In this work we use the 'standard' configuration with inner and outer walls maintained at constant temperatures T a and T b respectively, with no internal heating, and include a rigid lid (i.e., a no-slip upper boundary condition). The model works on a (R, φ, z) grid of 24 × 64 × 24 nodes stretched in R and z (to resolve the boundary layers). The parameters defining the annulus and working fluid (Table 1) are identical to the 'main comparison' of Hignett et al. [14], but use a range of Ω and ∆T = T b − T a , instead of just one value. Outer wall temperatures 20 Eighty-two simulations were run to construct a numerical regime diagram.
Each run was placed in the regime diagram using the thermal Rossby and Taylor numbers: where g is the gravitational acceleration and ν is the kinematic viscosity.
These parameters have been used for some decades as the standard quantities for comparing annulus flow regimes, and so will allow comparison between our results and earlier work. Both Θ and T are dimensionless parameters representing the balance of forces in the annulus. The thermal Rossby number Θ is the ratio between buoyancy/inertial forces and Coriolis forces; at large Θ rotation may be ignored, and at low Θ the geostrophic approximation is justified. The Taylor number T is the ratio between inertial forces due to rotation, and viscous forces; it measures the ability of a fluid to resist shear stresses. Using the values in Table 1, The simulations were initialized by first running a reduced model, which integrates the axisymmetric form of the equations of motion over a vertical 2D slice (with the same parameters as the subsequent 3D run). This was run until a steady state was reached (around t = 10, 000s). The slice was then copied to each azimuthal coordinate, and a sinusoidal perturbation was introduced

ACCEPTED MANUSCRIPT
to the temperature field, defined by where X = 0.1 K. Note that this perturbation is applied only at φ = 0, so it is essentially a white noise perturbation independent of any particular azimuthal behaviour. Each simulation starts from a different axisymmetric state (defined by the 2D run), but is initialized using the same perturbation.
The perturbation is necessary to move the system away from an (unstable) axisymmetric state, as it would otherwise remain axisymmetric (N.B.: the effects of using finite precision arithmetic would eventually grow to appreciable levels, causing the symmetry to be broken, but this would take a long time).
Each simulation was run at constant ∆T and Ω until transient motion had

ACCEPTED MANUSCRIPT
distinct signature; its trajectory in phase space is given by where ∆t is the time step size.

Dynamical regimes
Six distinct regimes were observed: axisymmetric flow (AX), steady wave flow (S), amplitude vacillation (AV), modulated amplitude vacillation (MAV), structural vacillation (SV) and a previously unseen regime, period-doubled amplitude vacillation (AV-d). The boundary between steady wave flow and vacillating flow was defined using the vacillation index [3], defined as To identify runs close to this boundary, two subdivisions of the steady wave regime have been defined, and the conditions for classifying runs in these  Table 2. S++ effectively defines a flow right on the boundary between steady wave flow and vacillating flow. Table 2 Classification of regimes based on vacillation index. I V is the mean value of I V over a time series (the mean of a set of I V calculated using each pair of adjacent maxima and minima) and σ I V is its standard deviation.
The main features of Fig and 'strong' otherwise; these are waves whose amplitudes are nonzero despite having a thermal Rossby number above the Eady threshold for baroclinic instability. One run was defined as being 'weak' in this context. The amplitude of the wave in this case was highest at mid-R and mid-z, in contrast with [19] whose weak waves were confined to the vertical extremes of the domain.

Differences between simulation and laboratory
A number of differences were also found between Fig. 3 and laboratory results.
Only a few flows exhibited modulated amplitude vacillation -in earlier laboratory work [3,18] it was more common, although the fluid was of higher

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
Prandtl number. We define MAV here as a flow where the dominant wave amplitude oscillates at one frequency but the oscillation is modulated; the time series spectrum contains a strong peak at the vacillation frequency and at least one peak at lower frequency. The attractor traced in the derivative coordinate embedding has nonzero thickness. Steady regimes were found at low ∆T , and period doubled amplitude vacillation was seen at high ∆T (Section 4). None of the more exotic regimes seen by Früh & Read [18] in fluid of higher Prandtl number were found.
Early laboratory work [19] found that, when changing Ω at constant ∆T during a single experiment, in some situations the regime transitions occur at different rotation rates when increasing Ω, compared with decreasing Ω (the so-called 'hysteresis' effect). In our work each run used constant Ω, so it is not yet possible to determine whether this effect also occurs in numerical simulations of the annulus.
Hignett [2] used the same annulus dimensions and fluid as these simulations, so our regime diagram is the equivalent numerical version of that laboratory study. In addition to the similarities in Section 3.1 above, in both laboratory and numerical regime diagrams a similar sequence of behaviour was ob-

ACCEPTED MANUSCRIPT
ratory setup [14]. Regime transitions are less clear at higher rotation rates, so it is impossible to say whether this shift is constant over the whole of (T , Θ) space. different from the amplitude vacillation in Fig. 2b. This difference is reinforced by comparing the frequency spectra of their amplitude time series: in Fig. 7b we see subharmonics below the main frequency, but these are not present in This defines a trajectory in phase space of length I = N − τ (D E − 1). A number of conditions restrict the governing parameters [21]: small sampling time t s = t A 2 − t A 1 ensures a smooth approximation to the true trajectory, the delay time τ should be long enough to ensure independence between coordinates, and the embedding dimension D E must be large enough to ensure a one-toone map from time series to delay reconstruction.

A sequence of bifurcations similar to the logistic map
In this work the sampling time was set such that 1.0 ≤ t s / s ≤ 2.5 (over- sampled, but at no extra computational cost). The embedding dimension was estimated using the method of false nearest neighbours [22], implemented using the Nonlinear Time Series Analysis (TISEAN) package [23,24]; D E = 4 was found to be suitable. The delay time was found by computing the first minimum of the time delayed mutual information [25], yielding values of 40 -120 t s (40 -300 s). The presence of these different periods indicates some kind of bifurcation sequence. To quantify this, for each run the set of temperature amplitude max- ima was identified in the T M =2 (t) time series, and these temperatures are plotted in Fig. 9 against the corresponding Taylor number for each run. The result is similar to behaviour produced, for example, by bifurcations of the logistic map [5][6][7]. At low Taylor numbers the behaviour is periodic and bifurcates to 2AV-d2 with increasing Taylor number. Further bifurcations occur before A (2AV-dh), whose values are spread over the temperature range. Two periodic windows (2AV-d3) separate A from B (2AV-dh). By C (2AV-d3) periodic behaviour has returned. Around T ≈ 4 × 10 6 an event occurs with results characteristic of a boundary crisis [28]. This causes the temperature range to tighten by a factor of 2-3, and periodicity is lost by D (2AV-d1).

ACCEPTED MANUSCRIPT
The 2AV-d3 cases may imply that chaotic dynamics exist at certain intervals of Taylor number along the sequence [29]. It is also notable that the run to the left of A is 2AV-d3: this indicates that within the range 3.388 × 10 6 ≤ T ≤ 3.508 × 10 6 there should be a chaotic window, as the first period-3 window

ACCEPTED MANUSCRIPT
for bifurcations of the logistic and similar maps occurs immediately after the first chaotic window.
We now consider whether a map similar to the logistic map exists to define the sequence of time series maxima in our simulations. In a continuous dynamical system dX/dt = F (X), the first return map is constructed from a Poincaré section. For instance, this section can be an (n − 1)-dimensional subspace defined where x n is constant: the first return map is then Y i+1 vs Y i . In this context, we define the section using x n ≡ dT M =2 /dt = 0. The first return map is a suitable fit, and the best fits are superimposed in all four cases. For 2AV-dh these fits are particularly good, indicating an underlying structure qualitatively similar to the logistic map.

A superposition of waves?
It is possible that this behaviour could be the produced by superimposing waves of incommensurate frequency, manifested in the dominant wavenumber at mid-height/mid-radius. We now argue that this is not the case, by showing that the oscillation of interest is present throughout the fluid, and is also present at other wavenumbers. To show that the oscillations are ubiquitous in space, the same spatial Fourier analysis was done at other points in the it would appear that this feature is typical of amplitude vacillating waves in the annulus.

ACCEPTED MANUSCRIPT
The ubiquity of the oscillations as a nonlinear wave-zonal flow interaction can be confirmed by analyzing an equivalent time series for the Nusselt number at one of the sidewalls (Fig. 10c, for R = a, the inner sidewall). The Nusselt number is a dimensionless parameter that quantifies the importance of convective over conductive heat transfer, and is a single value summed over the whole sidewall. The evolution of the Nusselt number in Fig. 10c  We now verify the connection with the logistic map using two further tests.  Figure   12 shows the first, second and third return maps generated from these data (on the left), and on the right the corresponding return maps of the logistic map using r = 3.8 (run A) and r = 3.87 (run B). For all three return maps, the similarities between the data and the logistic map are clear.

Map predictions
The ultimate test of whether the map is useful or not is its ability to predict the next maximum in the series. For each maximum in the A and B time

ACCEPTED MANUSCRIPT
series, the first return map fitted using Equation 1 was used to predict the next maximum in the sequence. The absolute error in these predictions was fairly uniformly spread between 0 K and 0.015 K with mean 0.0064 K for run A and between 0 K and 0.01 K with mean 0.0062 K for run B. Taking the size of the attractor to be the range in temperature covered by the maxima, this corresponds to an absolute mean error of 7% of the attractor size for run A, and 6% for run B. These values again confirm that the quadratic map is suitable to describe the sequence of temperature maxima. Finally, we verify that the 2AV-d flows in the bifurcation sequence are chaotic by computing estimates of their maximum Lyapunov exponents. If the dynamics follow the logistic map structure exactly, regimes 2AV-d1 through 2AV-d3

ACCEPTED MANUSCRIPT
should not be chaotic and 2AV-dh should be chaotic. Deterministic systems with a single positive Lyapunov exponent λ max > 0 that display aperiodic long-term behaviour satisfy the conditions for chaos [30, p.323]. Wolf et al. [31] define the ith Lyapunov exponent as In each case the exponent was measured by running the algorithm for embedding dimensions 3 ≤ D E ≤ 7; the maximum and minimum length scales were represented by 20% of the attractor diameter [peak to peak difference in T M =2 (t)] and 0.1% of the attractor diameter respectively, and a propagation time of 25% of the vacillation period was used. The exponent is the value that λ max (t) asymptotically tends towards, where t is the length of the time series used to define the attractor (as t → ∞ the whole attractor becomes available,

ACCEPTED MANUSCRIPT
hence the asymptotic value of λ max (t) is appropriate). If the plot did not contain such a region for this range in D E , the calculation was repeated with a longer delay time. The exponent was then taken to be the mean value over the dimensions and ranges in time where this was satisfied (the plot for 2048 points in Fig. 10 of Wolf et al. [31] gives an example of an appropriately flat region).
The exponent was calculated for each of the 28 runs, and the results are shown in Fig. 13, plotted against Taylor number on the same horizontal scale as Fig. 9.
The values are close to zero at low Taylor number, where period-1 and 2AV-d2 behaviour is observed. The rest of the plot is split into two regions of high λ max near A and D (where the aperiodic behaviour is observed, now confirmed as chaotic by these exponents), and a region of lower λ max near C, in the 2AV-d3 window. By comparison, λ max calculated for a selection of AV and MAV runs in Fig. 3 is in the range (−1.2 × 10 −4 ± 3 × 10 −5 ) ≤ λ max / bits s −1 ≤ (6.3 × 10 −4 ± 10 −4 ), and for the SV runs it is (1.3 × 10 −3 ± 7 × 10 −5 ) ≤ λ max / bits s −1 ≤ (5.5×10 −3 ±2×10 −4 ). It should also be noted that the runs in the periodic windows of the bifurcation sequence have Lyapunov exponents an order of magnitude larger than the exponents for the AV and MAV runs.
We can connect these results with the map predictions in Section 4.3. The λ max calculated here are of order 0.1 -1.0 bits / orbit, where the orbital period is the vacillation period of the dominant wave mode (e.g. Fig. 2b). The error-doubling time is therefore a few orbits, so we would not expect even a near-perfect map to predict successive maxima with perfect precision. In this light, a mean error of only ∼6% in the predictions is encouraging.
Behaviour resembling the logistic map has not been reported previously in the thermally-driven baroclinic rotating annulus, but there is a precedent elsewhere in the study of rotating fluids. Figures 6-8  Most of the literature on period-doubling in laboratory fluid systems is in the area of Rayleigh-Bénard convection [8][9][10][11]. Period-doubling is presented as one of a number of possible routes to chaos or turbulent convection, and is identified primarily through spectral analysis of time series. Similar analysis is possible here, but as the fine detail of parts of the sequence was not explored (particularly the bifurcations before the first chaotic window) this analysis is less conclusive. The appearance of subharmonics can be seen between Figs 2c and 7b: f /3 and f /1.5 are present in 7b but not in 2c. Evidence is limited elsewhere, however. Furthermore, the lack of runs in the region of the sequence where the initial bifurcations take place means that values of the Feigenbaum constants [6] cannot be calculated for our sequence. This is unfortunate, as they would allow further verification of the dynamics as well as comparison with the Rayleigh-Bénard system, where the data were accurate enough for these constants to be calculated.
It may be possible to observe period-doubled amplitude vacillation in the A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT
laboratory in real annulus experiments. At the most crowded point along our sequence, the points differ in rotation rate by ∆Ω ≈ 0.002 rad s −1 , and in temperature difference by ∆(∆T ) ≈ 0.05 K. By comparison, laboratory rotation rates are typically stable to ±0.0004 rad s −1 , and temperature differences to ±0.02 K [34]. It may also be easier to derive the Feigenbaum constants using laboratory results, if similar methods to the Rayleigh-Bénard work are used.
In conclusion, the evidence from Figs 5-8, 9, and 12 strongly supports an underlying structure similar to the logistic map. It adds to the already rich range of behaviour observed in the rotating annulus over the years. A better quantitative analysis of this regime would require more simulations run for longer, but due to the time constraints imposed by the simulation this is computationally very expensive. Finally, the bifurcations we have described here are much simpler than the regime transitions that occur in the dynamical equations themselves; we have found evidence for an underlying simplicity in these complicated flows, which is an important step towards understanding them more generally.
RMBY acknowledges support via NERC Studentship NER/S/A/2005/13667, and thanks A. A. Castrejón-Pita for a useful suggestion during revision. We also wish to thank two anonymous reviewers for their comments.