ABSTRACT
Title of dissertation: PROBLEMS IN SPATIOTEMPORAL CHAOS
Matthew Tyler Cornick
Doctor of Philosophy, 2007
Dissertation directed by: Edward Ott
Department of Physics
In this thesis we consider two problem areas involving spatiotemporally chaotic
systems. In Part I we investigate data assimilation techniques applicable to large
systems. Data assimilation refers to the process of estimating a system?s state from a
time series of measurements (which may be noisy or incomplete) in conjunction with
a model for the system?s time evolution. However, for practical reasons, the high di-
mensionality of large spatiotemporally chaotic systems prevents the use of classical
data assimilation techniques such as the Kalman fllter. Here, a recently developed
data assimilation method, the local ensemble transform Kalman Filter (LETKF),
designed to circumvent this di?culty is applied to Rayleigh-B?enard convection, a
prototypical spatiotemporally chaotic laboratory system. Using this technique we
are able to extract the full temperature and velocity flelds from a time series of
shadowgraphs from a Rayleigh-B?enard convection experiment. The process of es-
timating  uid parameters is also investigated. The presented results suggest the
potential usefulness of the LETKF technique to a broad class of laboratory experi-
ments in which there is spatiotemporally chaotic behavior.
In Part II we study magnetic dynamo action in rotating electrically conduct-
ing  uids. In particular, we study how rotation efiects the process of magnetic fleld
growth (the dynamo efiect) for a externally forced turbulent  uid. We solve the
kinematic magnetohydrodynamic (MHD) equations with the addition of a Corio-
lis force in a periodic domain. Our results suggest that rotation is desirable for
producing dynamo  ows.
PROBLEMS IN SPATIOTEMPORAL CHAOS
by
Matthew Tyler Cornick
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulflllment
of the requirements for the degree of
Doctor of Philosophy
2007
Advisory Committee:
Professor Edward Ott, Chair/Advisor
Professor Brian Hunt
Dr. Istvan Szunyogh
Professor Thomas Antonsen
Professor Daniel Lathrop
c Copyright by
Matthew Cornick
2007
Preface
This thesis is organized into two parts. Both are concerned with the study of
spatiotemporal chaos, a form of chaos characterized by disorder in both space and
time. Spatiotemporal chaos often arises in spatially extended systems which have a
flnite correlation length; in particular, those which have a correlation length much
smaller than the full system size. These complex systems are high dimensional; they
have many (often hundreds or thousands) of dynamical degrees of freedom. Exam-
ples of spatiotemporal chaos have been found in optics [1], chemical and biological
media [2], and hydrodynamics [3, 4], including geophysical  ows in the ocean and
atmosphere.
Part I is concerned with the problem of state and parameter estimation in
spatiotemporally chaotic systems. In particular we study the problem in relation to
a common experimental system, Rayleigh-B?enard convection. In Part II we report
on results from magnetohydrodynamic simulations of rotating turbulence.
Often spatiotemporal chaos is distinguished from turbulence [5, 6], and thus
many geophysical  ows, such as the atmosphere (which is highly turbulent) and
magnetohydrodynamic turbulence, would not be classifled as spatiotemporal chaos.
For our purposes this distinction is unimportant, but it is instructive to discuss the
difierences. Consider the various length scales in a system. The size of the full system
is denoted D while the correlation length is denoted ?c. In addition, the length scale
at which energy is introduced into the system is lE, while the length scale at which
energy is dissipated (usually the smallest scale in the system) is lD. When lD ?lE
ii
the system is said to be turbulent. While for lD ? lE, the system may exhibit low
dimensional temporal chaos (?c ? D) or spatiotemporal chaos (?c ? D). Despite
this distinction, both turbulence and spatiotemporal chaos are high dimensional,
spatially disordered states. Part II is concerned with a turbulent type of chaos,
while Part I is concerned mainly with Rayleigh-B?enard convection, an example of
non-turbulent spatiotemporal chaos. However, many of the results and techniques
described in Part I can (and have) been applied to turbulent systems (in particular,
for the purpose of weather forecasting).
iii
Dedication
For my wife, and Mr. Magney.
iv
Acknowledgments
Over the many years of my education there have been several people who
have have in uenced me to take the path towards a physics Ph.D. Two high school
teachers in particular, Mr. Magney and Mr. Rischling, were vital in my decision to
study Physics from the beginning. Their constant encouragement was a powerful
motivator. As an undergraduate, my advisor Dr. Field had a tremendous efiect on
my decision to study physics in graduate school. He constantly took interest in my
growth and understanding.
I owe my gratitude to the University of Maryland. I would especially like to
thank my advisor Edward Ott who provided his wisdom and experience throughout
my years at UMD. I would also like to thank him for his text on Dynamical Systems,
which, as a undergraduate, captured my interest and greatly in uenced my choice
of graduate school. Thanks to Brian Hunt, who gave me valuable advice on many
occasions. I am also tremendously grateful for the direction and advice provided by
Tom Antonsen and Dan Lathrop, for Mike Schatz who always made me feel welcome
during visits to Georgia Tech, and for Alex Dragt who gave me my flrst opportunity
at graduate research. Thanks to Nick Mecholsky for always having an interesting
problem to solve, and to Young-Noh Yoon, who constantly reminded me that even
the most basic things can be questioned. I am grateful to Laurette Tuckerman for
providing her Fortran code for simulating the Boussinesq equations.
Of course, without the encouragement of my parents, none of this would have
been possible. My father?s passion for science and nature drove me to similar inter-
v
ests. Lastly, thanks to my wife, who was with me at every step, and who reminds
me that the most important things cannot be found in a text book.
vi
Table of Contents
List of Tables ix
List of Figures x
List of Abbreviations xiii
I PARAMETER AND STATE ESTIMATION IN EXPERIMENTS EXHIBIT-
ING SPATIOTEMPORAL CHAOS 1
1 Introduction 2
1.1 Rayleigh-B?enard Convection . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.1 The Shadowgraph Method . . . . . . . . . . . . . . . . . . . . 12
1.2 Data Assimilation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3 The Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.3.1 Bayesian Estimation . . . . . . . . . . . . . . . . . . . . . . . 16
1.3.2 Kalman Filter Equations . . . . . . . . . . . . . . . . . . . . . 18
1.3.3 The Extended Kalman Filter . . . . . . . . . . . . . . . . . . 19
1.4 Ensemble Kalman Filters . . . . . . . . . . . . . . . . . . . . . . . . 20
1.4.1 Parameter Estimation . . . . . . . . . . . . . . . . . . . . . . 25
1.4.2 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2 Application of Data Assimilation to Rayleigh-B?enard Convection 27
2.1 The Local Ensemble Transform Kalman Filter . . . . . . . . . . . . . 29
2.2 Direct Insertion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3 Results 36
3.1 Perfect Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.1.1 Performance with Noise/Sparseness . . . . . . . . . . . . . . . 38
3.1.2 Parameter Estimation . . . . . . . . . . . . . . . . . . . . . . 43
3.2 Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.2.1 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . 47
3.2.2 Analysis of Measurement Noise . . . . . . . . . . . . . . . . . 48
3.2.3 Parameter Estimation and Performance with Sparseness . . . 51
3.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
II THE EFFECT OF ROTATION ON DYNAMO ACTION IN MAGNETO-
HYDRODYNAMIC TURBULENCE 63
4 Introduction to Magnetohydrodynamics 64
4.1 The Kinematic Dynamo . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.2 Rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
vii
5 Numerical Approach 73
5.1 Spectral Representation . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.2 Time Stepping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.3 Parameters and Resolution Requirements . . . . . . . . . . . . . . . . 78
6 Simulation Results 80
6.1 Flow Characterization . . . . . . . . . . . . . . . . . . . . . . . . . . 82
6.1.1 Flow Structure . . . . . . . . . . . . . . . . . . . . . . . . . . 90
6.1.2 The Energy Spectrum . . . . . . . . . . . . . . . . . . . . . . 94
6.2 Kinematic Dynamo Characterization . . . . . . . . . . . . . . . . . . 97
6.3 Connection to an Envisioned Experimental Situation . . . . . . . . . 103
6.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
6.4.1 Future Prospects . . . . . . . . . . . . . . . . . . . . . . . . . 110
A Geometric Optics Derivation of the Shadowgraph Light Intensity 111
B Short-Time Nonlinear Evolution of an Initially Gaussian PDF 114
Bibliography 121
viii
List of Tables
6.1 The average Urms, R, Ro, and Ek for various ?, ??1 = 36. . . . . . . 81
6.2 The average Urms, R, Ro, and Ek for various ?, ??1 = 18. . . . . . . 81
6.3 The average Urms, R, Ro, and Ek for various ?, ??1 = 6:3. . . . . . . 82
ix
List of Figures
1.1 Images from simulated Rayleigh-B?enard convection, demonstrating
the time evolution of the system state. . . . . . . . . . . . . . . . . . 7
1.2 The z-dependence of the temperature fleld. . . . . . . . . . . . . . . . 7
1.3 Sequence of images of the temperature error. . . . . . . . . . . . . . . 11
2.1 Local regions on the cylindrical mesh. . . . . . . . . . . . . . . . . . . 30
3.1 Temperature error of perfect model test under ideal conditions. . . . . 39
3.2 Mean Flow error of perfect model test under ideal conditions. . . . . 40
3.3 Quality measures (Emin , Emin?u , and the predictability time ?) as the
density of observations is reduced . . . . . . . . . . . . . . . . . . . . 41
3.4 Quality measures (Emin , Emin?u , and the predictability time ?) as mea-
surement noise is increased. . . . . . . . . . . . . . . . . . . . . . . . 42
3.5 Demonstration of the convergence of the parameters in perfect model
tests. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.6 Parameter estimation of R in a small aspect ratio test, comparison
between the LETKF and the EKF. . . . . . . . . . . . . . . . . . . . 46
3.7 Parameter estimation of Pr in a small aspect ratio test, comparison
between the LETKF and the EKF. . . . . . . . . . . . . . . . . . . . 46
3.8 Probability distribution of shadowgraph pixel noise. . . . . . . . . . . 49
3.9 Correlation of experiment shadowgraph pixel noise to nearby pixels. . 50
3.10 Correlation C(m;n) of measurement noise. . . . . . . . . . . . . . . . 51
3.11 Image of mean  ow vectors overlaying a experimentally obtained
shadowgraph. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.12 An estimate of the  uid state showing a shadowgraph measurement,
the temperature proflle, the modeled shadowgraph, and the inferred
vorticity potential. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.13 Forecast error EI for DI and LETKF methods are shown for high and
low measurement densities. . . . . . . . . . . . . . . . . . . . . . . . . 55
x
3.14 Forecast error EI for DI and LETKF methods in perfect model (PM)
tests and when using experimental data (E). . . . . . . . . . . . . . . 56
3.15 E dimension decreasing as the LETKF converges on perfect model
data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.16 E dimension decreasing as the LETKF converges on experimental data. 60
3.17 Sequence of images showing the convergence process when a large
data void is present in shadowgraph measurements. . . . . . . . . . . 61
6.1 Time series of the force magnitudes: inertia, pressure, viscous, and
Coriolis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
6.2 Time series of the velocity magnitude, helicity, and power. . . . . . . 84
6.3 Average velocity versus rotation rate. . . . . . . . . . . . . . . . . . . 85
6.4 Average inertia versus rotation rate. . . . . . . . . . . . . . . . . . . . 86
6.5 Average pressure force versus rotation rate. . . . . . . . . . . . . . . . 86
6.6 Average power versus rotation rate. . . . . . . . . . . . . . . . . . . . 87
6.7 Average helicity versus rotation rate. . . . . . . . . . . . . . . . . . . 87
6.8 Average viscous force versus rotation rate. . . . . . . . . . . . . . . . 88
6.9 Average Coriolis force versus rotation rate. . . . . . . . . . . . . . . . 88
6.10 Power spectrum showing presence of inertial waves. . . . . . . . . . . 90
6.11 Images of the pressure on the cube surface. . . . . . . . . . . . . . . . 91
6.12 Images of the velocity magnitude on the cube surface. . . . . . . . . . 92
6.13 Images of the Coriolis force magnitude on the cube surface. . . . . . . 93
6.14 Images of the velocity magnitude on the cube surface for ??1 = 36
(N=64). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
6.15 High resolution kinetic energy spectrum, validating the code. . . . . . 95
6.16 Compensated energy spectrum for various rotation rates. . . . . . . . 96
6.17 The magnetic energy growth rate  ? as a function of ? for several
initial conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
xi
6.18 Growth rate versus magnetic Reynolds number. . . . . . . . . . . . . 99
6.19 The critical ? versus ?. . . . . . . . . . . . . . . . . . . . . . . . . . . 100
6.20 The critical Rm versus ?. . . . . . . . . . . . . . . . . . . . . . . . . . 101
6.21 Images of the magnetic fleld strength on the cube surface. . . . . . . 102
6.22 Surface seperating dynamo behavior in the simulations. . . . . . . . . 105
6.23 Curve separating dynamo behavior from non dynamo behavior in the
plane (P;C) for Prm = 0:5. . . . . . . . . . . . . . . . . . . . . . . . 106
6.24 The critical power curves P(??)c , P(?)c , and PTc for Prm = 0:2. . . . . 108
6.25 Critical power curves P(+)c and P(?)c versus C for several values of Prm.109
xii
List of Abbreviations
KF Kalman fllter
EKF extended Kalman fllter
EnKF ensemble Kalman fllter
LEKF local ensemble Kalman fllter
LETKF local ensemble transform Kalman fllter
RMS root mean square
MHD magnetohydrodynamics
FFT fast Fourier transform
CFL Courant-Friedrichs-Levy
xiii
Part I
PARAMETER AND STATE ESTIMATION IN EXPERIMENTS EXHIBITING
SPATIOTEMPORAL CHAOS
1
Chapter 1
Introduction
Estimation of the state of an evolving dynamical system from measurements is
often a prerequisite for prediction and control of the system. However, obtaining the
system state is a common experimental di?culty for many spatiotemporally chaotic
systems, where available measurements may be incomplete and noisy. When an
approximate model for the system is available, it can be used in conjunction with
incoming measurements to estimate the evolving system state, a process referred to
as ?data assimilation?.
Most data assimilation algorithms are iterative, cycling between a predict and
update step once every time interval ?t. In the update step, current measurements
are used to update (or correct) the prediction. The prediction step then propagates
the updated state, via the model, to the next measurement time (i.e., it is a short
term forecast). The aim of this process is to synchronize the system and its model
by coupling them via the measurements.
The Kalman fllter [7], described in Section 1.3, optimally solves the data as-
similation problem for systems with linear dynamics. Several methods extending the
Kalman fllter methodology to nonlinear systems have been proposed, including the
extended Kalman fllter (EKF) [8] (Section 1.3), and the class of ensemble Kalman fll-
ters (EnKF) [9] (Section 1.4). Straightforward application of these methods to large
2
spatiotemporally chaotic systems is often completely infeasible. In particular, the
EKF requires inversion ofN?N matrices, whereN is the number of model variables.
This can be quite prohibitive; for a discretized partial difierential model evolving
M scalar spatial flelds in time, the number of model variables is M multiplied by
the number of grid points (which can number in the millions). In addition, the en-
semble (EnKF) methods require prohibitively large ensemble sizes when applied to
high-dimensional systems. Despite these di?culties, recent developments [10, 11, 12]
from the fleld of numerical weather prediction [13, 14, 15, 16, 17, 18] suggest the
possibility of achieving good accuracy (as in a Kalman fllter), but in a way that is
computationally feasible for large nonlinear systems.
In this thesis we demonstrate the e?cacy of a new method, the local ensemble
transform Kalman fllter (LETKF) [12] (Chapter 2). Although originally motivated
by application to weather prediction, the LETKF is potentially broadly applicable to
any large spatiotemporally chaotic system. In particular, we investigate Rayleigh-
B?enard convection (Section 1.1). Flows such as spiral defect chaos [19, 3] in the
Rayleigh-B?enard problem are, perhaps, the best studied experimental examples of
spatiotemporal chaos; nevertheless, many general aspects of spiral defect chaos re-
main poorly understood.
The LETKF is motivated by the observation that spatiotemporally chaotic
systems exhibit a flnite correlation length much smaller than the system size. Fre-
quently, regions much smaller than the system size can be described by relatively
few degrees of freedom. With this in mind, the LETKF employs many independent
data assimilations in a set of overlapping local regions, each with a characteristic
3
length on the order of the correlation length. Because these regions are relatively
small, their individual computations are not prohibitive. Furthermore, by use of a
simple example [10, 11] it was indicated that, by exploiting localization in this way,
state estimates with accuracies virtually the same as those for a classical Kalman
fllter technique (thus presumably of near optimal accuracy) can be achieved.
What follows is an introduction to Rayleigh-B?enard convection in Section 1.1
and an introduction to the classical methods of data assimilation in Sections 1.2, 1.3,
and 1.4. Details of the LETKF algorithm, as well as our application of the LETKF
to Rayleigh-B?enard convection is described in Chapter 2. Tests of the accuracy of the
LETKF are presented in Sections 3.1 and 3.2, in which we investigate performance
with extremely sparse/noisy measurements and test extensions of the LETKF for
estimating model parameters.
1.1 Rayleigh-B?enard Convection
In Rayleigh-B?enard convection, a horizontal  uid layer of thickness d is con-
flned between a heated lower plate and a cooled upper plate. For a temperature
difierence between the plates ?T that is su?ciently small, the  uid is at rest, and
heat is transported by conduction. In this state, the temperature rises linearly from
the lower boundary with temperature TH to the upper boundary with temperature
TC = TH ??T. The onset of  uid motion occurs when buoyancy overcomes viscous
dissipation and thermal difiusion as ?T is raised above a critical value ?Tc.
Initial analytic work on Rayleigh-B?enard convection showed that, in a space
4
of inflnite horizontal extent, the system has a stable state consisting of parallel
convection rolls. However, Rayleigh-B?enard convection was recently shown [3] to
support a type of spatiotemporal chaos which has been named spiral defect chaos
due to the abundance of spiral structures and roll defects present in the evolving
state. This state has been extensively studied theoretically, numerically, and exper-
imentally [20, 21, 22, 23, 24, 25, 6, 3, 26]. For a recent review of Rayleigh-B?enard
convection see [19].
Rayleigh-B?enard convection is typically modeled using the Boussinesq equa-
tions [27], which are commonly nondimensionalized with temperature scaled by ?T,
length scaled by d, and time scaled by the vertical difiusion time tv = d2=?, where
? is the thermal difiusivity. This system of units is used throughout Part I of this
thesis. The temperature fleld is denoted T, while the temperature deviation from
the conducting static solution is denoted  ; T(x;y;z) = ?T (x;y;z) +TH ? ?T zd.
We solve the Boussinesq equations in the disk shaped region x2 +y2 ? ?2, jzj? 12,
with Dirichlet boundary conditions u = 0,  = 0 on all walls. ? is the radius of the
disk in units of d and is often referred to as the aspect ratio. The  boundary con-
dition represents the situation which would occur for a conducting wall (when the
wall?s thermal conductivity much higher than the  uid?s thermal conductivity). In
terms the  uid?s velocity u, temperature deviation  , and pressure p, the Boussinesq
equations take the form
 @
@t + u?r
?
u = ?rp+Prr2u +PrR ^z ;
 @
@t + u?r
?
 = r2 + u? ^z ; (1.1)
5
r?u = 0 :
These equations have two dimensionless parameters, the Rayleigh number R and
the Prandtl number Pr,
R = gfid
3?T
?? ;Pr =
?
?: (1.2)
Here fi is the thermal expansion coe?cient, ? is the kinematic viscosity, and g is
gravitational acceleration (in the ?^z direction). The critical Rayleigh number for
convective onset is Rc ? 1707. The reduced Rayleigh number
? = R?RcR
c
= ?T ??Tc?T
c
(1.3)
measures the amount above onset. Fluid convection arises when ?> 0.
We have investigated the parameter region near ? = 1; Pr = 1. At these
values of ? and Pr, the spatiotemporally chaotic state known as spiral defect chaos
can arise [3, 19]; however, in our studies using ? ? 20, the region is too small to
support the large spirals typically seen in spiral defect chaos. Nevertheless, the
convective  ows in our studies exhibit complex behavior in both space and time.
See Fig. 1.1 for an example of the spatial structure of the evolving state (the images
show simulated shadowgraphs, an imaging technique described in Section 1.1.1).
6
Figure 1.1: Images from simulated Rayleigh-B?enard convection (? = 1, ? = 20),
demonstrating the time evolution of a typical system state. White represents cold
descending  uid, while black represents warm rising  uid.
The z-dependence of the  (x;y;z) fleld from a typical state on the system
attractor exhibits asymmetry about z = 0, as shown in Fig. 1.2.
Figure 1.2: The z-dependence of the temperature deviation is shown in this cross-
section (the x = 0 plane). Red indicates that the temperature is higher than the
conducting proflle, whereas blue indicates that it is lower.
Our technique for integrating equations (1.1) is the pseudospectral method
7
described in [28]. Brie y, it uses a backward Euler time step for the linear terms,
and a second order Adams-Bashforth method for the nonlinear terms. All flelds are
expressed in terms of their Fourier components in the periodic azimuthal ` direction,
and in terms of Chebyshev polynomials in the bounded directions ?1=2 <z < 1=2
and 0 < r < ?. Thus the physical space representation exists on a mesh of grid
points (rm; `n; zl) which are evenly spaced in ` and located at the Chebyshev
points in z and r; `n = 2?nN
`
; rm = ? cos( ?(m?1)2(Nr?1)); zl = 12 cos(?(l?1)Nz?1 ). Here, the
number of grid points in the r, z, and ` direction are Nr, Nz, and N` respectively.
The dynamics of this system are driven by two instabilities. The skewed
varicose instability, for large wave numbers (small wavelengths) and the Eckhaus
instability (dilation instability), associated with small wave numbers (large wave-
lengths). Between these extremes lie stable wave numbers inside the Busse balloon,
a compact region of stable wavelengths in the space (?; k). The motion of defects
in the roll planform are ultimately a consequence of these instabilities.
Small structures (defects) rapidly evolve on a time scale of approximately one
vertical difiusion time. Large structures evolve on a time scale of one horizontal
difiusion time th = ?2tv. The correlation length ? decreases with increasing ? ap-
proximately according to a power law, ? ???1=2 [19]. For the investigated parameter
region (Pr; ? ? 1) the correlation length of this system is much smaller than the
diameter of the domain; approximately 2d to 3d.
Other than the Boussinesq equations, several models for Rayleigh-B?enard con-
vection have been proposed and studied [6]. In particular, the modifled Swift-
Hohenberg equations [29, 30] reduce the problem to two dimensions; compared to
8
three dimensions for the full Boussinesq equations (and thus they may be simu-
lated more easily than the full Boussinesq equations). These equations are used to
model many spatiotemporally chaotic phenomena, including spatiotemporal chaos
in lasers [31]. They are
@?
@t +gmu?r? =
h
?0 ?(1 +r2)2
i
???3; (1.4)
@v
@t =
?@?
@y
@f
@x ?
@?
@x
@f
@y
?
+Pr?r2 ?c2?v: (1.5)
The ? fleld is similar to the midplane temperature fleld of the Boussinesq equations,
the v = r2? fleld is vertical vorticity of the  uid  ow u = @?@y ^x? @?@x ^y, and f = r2?.
The parameter ?0 is related to ? via the relation ?0 ? 2:78?. These equations, when
simulated with boundary conditions ? = v = ? = f = 0 on the boundary of a
circular region, result in patterns very similar to those found in experiments on
Rayleigh-B?enard convection. In particular, large rotating spirals and defects are
present. Although these equations produce qualitatively similar solutions to the
Boussinesq equations, they exhibit spiral defect chaos only as a transient behavior.
Thus, in contrast to the full Boussinesq equations, spiral defect chaos does not occur
on the Swift-Hohenberg attractor. Because of this drawback we decided to use the
full Boussinesq equations as our model of Rayleigh-B?enard convection, despite the
computational complexity of simulating a full three dimensional domain.
Spiral defect chaos in Rayleigh-B?enard convection exhibits sensitivity to initial
conditions. It is instructive to demonstrate the error growth of an initial perturba-
tion. To do so, we simulate two  uid states, one termed the reference, and the other
which is identical to the reference state except for a small local perturbation in the
9
initial condition. The state of the reference fleld is denoted Tr(x;y;z);ur(x;y;z),
while the perturbed state is denoted T(x;y;z);u(x;y;z). Fig. 1.3 shows a sequence
of images demonstrating the nonlinear growth of T(x;y;0)?Tr(x;y;0) (in the mid-
plane z = 0). The error initially decays, and the two states remain close for several
tv. When the states diverge it is usually at the location of a defect. For example,
in the reference state a roll may be pinched ofi by the action of the skewed varicose
instability forming a defect while the second state may not undergo the pinch-ofi
event. This defect related error growth quickly spreads across the entire disk in
a time much less than the horizontal difiusion time. This highlights an important
property of spiral defect chaos; the dynamics are stable to small local perturbations
almost everywhere. It is only near some defects that the perturbation grows. This
was shown in [20]. The Lyapunov vector associated with the largest Lyapunov ex-
ponent is negative almost everywhere (except near a few defects, where it is large
and positive).
10
Figure 1.3: Sequence of images showing the growth of a small initial perturbation
located in the upper center of the disk at frame 0. Contours are the curves for
which Tr(x;y;0) = 0. Red indicates a slightly warmer midplane temperature than
the reference state, while blue indicates a slightly cooler temperature. The time
spanned between frame 0 and frame 6 is approximately 10tv.
11
1.1.1 The Shadowgraph Method
In experiments, Rayleigh-B?enard  ows are visualized using the shadowgraph
method [32]. The typical setup consists of a re ecting bottom plate and a trans-
parent top plate, conflning the (transparent)  uid. Light shines down through the
 uid, re ects from the bottom plate, travels up through the  uid and out of the disk
region. The exiting light is modifled by the index of refraction of the  uid (which has
a slight temperature dependence) and imaged via a CCD camera. With this method
the vertically integrated temperature fleld of the  uid is indirectly measured.
We connect the shadowgraph light intensity I(x;y) to the temperature fleld
in the  ow using the relation derived from geometric optics [32, 33]
I(x;y) = I?(x;y)1?ar2
?? (x;y)
: (1.6)
Here, r2? ?@2=@x2 +@2=@y2 is the horizontal Laplacian, and the temperature fleld
is vertically averaged: ? (x;y) ?R  (x;y;z)dz. I?(x;y) is the incident light intensity
and a = 2z1jdn=dTj, where n is the index of refraction of the  uid, z1 is the optical
path length from the midplane of the  uid layer to the image plane (in units of d),
and the temperature coe?cient of the index of refraction jdn=dTj is evaluated at
the average temperature of the  uid layer. Equation (1.6) is derived in Appendix A.
For geometric optics to be valid, the condition kar2?? k? 1 must hold.
This remote sensing method represents an incomplete observation of the full
 uid state. Thez-dependence of the temperature fleld is unmeasured, as is the entire
velocity fleld. The velocity fleld is largely composed of the circulating roll velocity.
However, there is a smaller dynamically important component to the velocity fleld
12
known as the mean  ow ?u(x;y) (which advects the roll pattern),
?u(x;y) ?
Z
u?(x;y;z) dz; (1.7)
where u = u? +uz ^z. This component of the velocity fleld plays a signiflcant role in
the dynamics [24] and is di?cult to measure in the experiments. In order to achieve
prediction and control of Rayleigh-B?enard convection, accurate estimates of the  uid
state are required, including the mean  ow. Thus we pursue the use of data assimi-
lation techniques for the purpose of estimating the  uid state (u(x;y;z);  (x;y;z))
from a sequence of shadowgraph images.
1.2 Data Assimilation
In the most general situation, we are given a time series of measurements from a
particular dynamical system, where the measurements are made at the times tj (j =
1;2:::). These measurements may not describe the system state ? exactly; they may
be incomplete, sparse, or noisy. Suppose that, at each time, a measurement may be
written as a collection of scalar quantities so that we can represent a measurement
by a s-component vector y, and the time series as the collection of the vectors
y1; y2; :::. Our goal, at each time tj, is to infer the N component system state ?j
from the time series y1; y2; :::.
We begin with the dynamical model,
?j+1 = G(?j); (1.8)
which describes the time evolution of the system from a time tj to tj+1. In this thesis
we shal not consider the situation where there is an added random component on
13
the right hand side of (1.8). We assume the system evolves in continuous time,
d?
dt = F(?); (1.9)
so that G is an integration of (1.9) from t = tj to tj+1. G may be exact, though
generally it describes only approximately the dynamics of the true system. This
discrepancy, when it exists, is known as model error. In addition, we assume that
measurements are ideally (in the absence of measurement noise) uniquely determined
by the system state, y = M(?). M is known as an observation operator, the
mapping of ? into the observation space; it need not be a one-to-one mapping.
When measurement noise is present, individual measurements take the form yj =
M(?j) + ?. Here, ? represents a random error, which is assumed to be normally
distributed with mean zero and (s?s) covariance matrix R. In addition to model
error, the observation operator M may be approximate. We refer to this form of
error as observation operator error.
There are many ways to view the above problem. In one sense it can be
thought of as an inversion problem: flnd ?1; ?2; ::: given y1; y2; :::. Even though
M may not be invertible, this is possible because all measurements can be used, in
conjunction with the model, to determine each ?j. In other words, each individual
yj may not contain enough information to determine the state ?j, but the entire
time series y1; y2; ::: contains quite a bit of information, potentially enough to
obtain any individual ?j.
In a difierent sense, we can think of the problem as flnding the trajectory in
state space satisfying (1.8) that, when projected into the observation space, best flts
14
the measured time series. In this sense we can construct a quadratic cost function
to be minimized in order to obtain a maximum likelihood estimate of the trajectory.
Since we have assumed measurement noise follows Gaussian statistics, the most
sensible cost at each time is given by c(?;y) = (M(?)?y)T R?1(M(?)?y). To see
why this form is used, consider the vector ? = S?1(M(?j)?yj) where S is real and
symmetric and solves S2 = R. We havec(?;y) = (M(?)?y)T S?1S?1(M(?)?y) =
?T? = j?j2 . Roughly speaking, the cost is then the square ?distance? between
? and y in units of the number of standard deviations they are separated in the
observation space.
The trajectory is uniquely determined from the state at ant particular time
tk: ?j = Gj?k(?k) 8j, where the notation Gn(?) means applying G to ? n times (G
is invertible so n may be negative). Using this notation, the cost of a trajectory is
given by
C(?k) =
X
j
c(Gj?k(?k);yj);
where the trajectory is constrained to pass through the point ?k at time tk.
Yet another perspective would be to look at the problem as a synchronization
problem, as alluded to in the introduction. If we can use the measurements to
synchronize the model with the true system, we will have obtained a trajectory
similar to the one sought above.
The last two perspectives, that of synchronization, or flnding the trajectory
with maximum likelihood, are instructive. They propose not only ways of thinking
about the problem, but ways of solving it as well, i.e. by synchronizing the model
15
to the true system, or by minimizing a cost function. When G and M are linear,
there is a solution which accomplishes all the above objectives, the Kalman fllter.
1.3 The Kalman Filter
Suppose that the model and observation operator are linear, G(?) = Q?, and
M(?) = H?, where H and Q are matrices of sizes?N andN?N, respectively. The
Kalman fllter is an recursive Bayesian approach to the problem of state estimation.
It is provably optimal given the statement of the problem (Gaussian measurement
noise, linear G and M). Before describing the Kalman fllter, we introduce the
concepts of Bayesian estimation.
1.3.1 Bayesian Estimation
Suppose we have a probability distribution function (PDF) p(?) representing
our knowledge of the system state. It is unimportant how exactly we came about
this information. When a new piece of information comes along, in the form of a
measurement y, we wish to update the PDF to re ect the new information. Denote
the conditional probability of measuring y, given that the state is ?, asp(yj?). Using
Bayes? rule, the new PDF of the state, re ecting the new information provided by
the measurement is then
p(?jy) = p(yj?)p(?)p(y) ; (1.10)
16
where p(y) = R p(yj?)p(?)d? is simply a normalizing factor. If we make the as-
sumption that p(?) is Gaussian, we have
p(?) ? exp
?
?12(? ? ??p)TP?1(? ? ??p)
?
; (1.11)
for some covariance matrix P and mean ??p. The assumption that measurement
noise is Gaussian leads to
p(yj?) ? exp
?
?12(y ?H?)T R?1(y ?H?)
?
: (1.12)
Using equation (1.10) we get the updated PDF
p(?jy) ? exp
?
?12 ?(y ?H?)T R?1(y ?H?) + (? ? ??p)TP?1(? ? ??p)?
?
: (1.13)
This PDF is Gaussian as well. Expressing p(?jy) in the form
p(?jy) ? exp
?
?12(? ? ??u)T U?1(y ? ??u)
?
; (1.14)
leads to the relation
(????u)T U?1(????u)+C? = (y?H?)T R?1(y?H?)+(????p)TP?1(????p): (1.15)
The constant C? may take any value since it leavs p(?jy) proportional to p(yj?)p(?)
(any normalization is flxed by the factor p(y)). In order for this relation to hold for
all ?, we must have
U = ?P?1 + HT R?1H??1; (1.16)
??u = U ?P?1??p + HT R?1y?: (1.17)
To summarize, beginning with an initial Gaussian PDF p(?) with mean ??p and
covariance P, we have included new information obtained from a measurement made
17
in the observation space with mean y and covariance R. The end result is the
updated PDF with mean ??u and covariance U, given by equations (1.16) and (1.17).
1.3.2 Kalman Filter Equations
With the tools of Bayesian estimation established, we proceed with the Kalman
fllter formulation. Note that, since G is linear, that a Gaussian PDF transforms,
under the action of G to another PDF with Gaussian form. The PDF with mean
??uj?1 and covariance Uj?1 at a time tj?1 transforms into the PDF at time tj with
mean and covariance
??pj = Q??uj?1; (1.18)
Pj = QUj?1QT: (1.19)
This PDF is known as the predicted PDF; it is then updated using (1.16) and (1.17),
forming the updated PDF with mean and covariance,
Uj = ?P?1j + HT R?1H??1; (1.20)
??uj = Uj ?P?1j ??pj + HT R?1yj?: (1.21)
These update equations are often rewritten in the form
Uj = ?I + PjHT R?1H??1 Pj; (1.22)
??uj = ??pj + UjHT R?1(yj ?H??pj): (1.23)
Equations (1.18), (1.19), (1.22), and (1.23) are the complete Kalman fllter
equations. In words, the Kalman fllter evolves a PDF of the system state forward
18
in time with successive updates whenever measurements are made. This process of
predict/update steps continues in an iterative fashion for all measurements in the
time series. At any given time, the PDF represents the accumulated information
provided by all previous assimilated measurements. The process typically begins
with a covariance P1 which is large enough that the initial PDF spans the entire
system attractor (or the range of possible system states).
1.3.3 The Extended Kalman Filter
In the case that G and M are not linear, we can still perform the Kalman fllter
steps by linearizing them about the current state estimate in equations (1.19), (1.22),
and (1.23). Although this allows one to proceed with the computation, it is no longer
provably optimal. However, if the nonlinearities are weak, or measurements are very
frequent, linearization may be a successful strategy. We have the linear maps
^G(?;??) = G(??) + DG(??)(? ???); (1.24)
^M(?;??) = M(??) + DM(??)(? ???); (1.25)
where DG(??) is the Jacobian of G evaluated at the point ??, and DM(??) is the
Jacobian of M evaluated at the point ??.
Replacing the Kalman fllter equations with
??pj = G(??uj?1); (1.26)
Pj = DG(??uj?1)Uj?1DG(??uj?1)T; (1.27)
Uj = ?I + PjDM(??pj)T R?1DM(??pj)??1 Pj; (1.28)
??uj = ??pj + UjDM(??pj)T R?1(yj ?M(??pj)); (1.29)
19
gives the extended Kalman fllter (EKF) equations.
Since the EKF does not treat the nonlinear evolution of the PDF, it is of
interest how quickly the exact PDF (if it were evolved with the full nonlinear model)
would deviate from a normal distribution. The growth in the skewness and kurtosis
of the distribution measure this tendency to deviate from normality, which gives a
sense for the accuracy of the approximation made in (1.27). Appendix A treats this
problem by considering the initial time evolution of the skewness and kurtosis of a
Gaussian PDF evolving in a generic nonlinear model. The growth of the standard
skewness and kurtosis measures is shown to be second order in time for short times.
One method for roughly compensating for the nonlinear evolution of the covariance
matrix is the method of multiplicative variance in ation, in which the predicted
covariance matrix is in ated by a factor ?2 > 1 [16, 12].
Although the EKF has solved (approximately) the problem of state estimation
for nonlinear systems, it still must compute the inverse of at least one N?N matrix.
When N numbers in the millions this is computationally infeasible. In addition,
consider that with 8 bytes of precision per  oating point number, simply storing P
when N = 106 would require 4 terabytes (even after taking advantage of symmetry).
These restrictions limit the applicability of the EKF.
1.4 Ensemble Kalman Filters
For a spatially extended system evolving according to a set of partial difier-
ential equations, the system state is usually represented as a collection of values
20
(for each variable of interest) on each grid point of a mesh. This collection of vari-
ables may number in the millions (N ? 106). Often, the system investigated has
many fewer dynamical degrees of freedom, in the sense that the dimension of the
attractor is much less than N. One could conceivably take advantage of this fact
by constructing a reduced rank representation of the covariance matrices U and P.
The class of ensemble Kalman fllters (EnKF) accomplish this through the
representation of the PDF, not by its mean and covariance, but by an ensemble of
system states. This ensemble gives a flnite sampling approximate representation of
the probability distribution function (PDF) of the system state. For an ensemble
?1; ?2; ::: ?k with k members, the PDF mean ?? is
?? = 1
k
kX
i=1
?i; (1.30)
while the covariance C is given by
C = 1k?1ZZT; (1.31)
where the matrix Z has columns
Z = ???1j??2j:::j??k?; (1.32)
and ??i are the ensemble perturbations
??i = ?i ? ??: (1.33)
The EnKF approach is to approximate the PDF as lying in the space spanned by the
k ensemble members [9, 34, 16, 14]. The number of ensemble members must then be
large enough to account for the many dynamical degrees of freedom in the system,
21
but it need not be so large as to span the entire N dimensional space of possible
states. The general procedure is to compute the updated ensemble f?u;1::: ?u;kg
from an update of the predicted ensemble f?p;1::: ?p;kg,
update step : f?p;1j ::: ?p;kj g+
fmeasurementsg!f?u;1j ::: ?u;kj g (1.34)
predict step : ?p;ij+1 = G(?u;ij ) i = 1:::k: (1.35)
This iterative procedure begins with an initial predicted ensemble f?p;10 ::: ?p;k0 g
consisting of states randomly sampled from the system attractor. The maximum
likelihood estimate of the system?s state after an update step is the center of the
updated ensemble given by (1.30). Here, the entire ensemble evolves nonlinearly,
not just the mean as in equation (1.26) of the EKF. Furthermore, in the space
spanned by the ensemble, only k?k matrices are required, rather than the N ?N
matrices of the EKF. Although the EnKF comes in many forms, they all approach
the problem in this way. Here we follow a particular type of EnKF known as an
ensemble square-root Kalman fllter [14].
The predicted ensemble ?p;i has mean ??p and covariance P, deflned as in (1.30)
and (1.31). From this point on, all quantities are assumed to have a subscript j (we
are performing the t = tj update step). We wish to flnd the updated ensemble ?u;i
with mean ??u and covariance U, also deflned as in (1.30) and (1.31). This update
step is performed in the space of k-component vectors ~? which transform to the full
state space via
? = ??p + Zp~?; (1.36)
22
where Zp is deflned as in (1.32) for the predicted ensemble ?p;i. In this space of ~?
vectors, the predicted ensemble has mean zero and a covariance which is simply a
multiple of the k?k identity matrix, ~P = (k?1)?1I.
Each predicted ensemble member is projected into the observation space by
forming the ensemble yp;i = M(?p;i), using the full (possibly nonlinear) form of
M. The observation mean ?yp is deflned similarly to (1.30). In addition, we deflne
the matrix Y p, containing the ensemble perturbations, Y p = ??yp;1j?yp;2j:::j?yp;k?,
where ?yp;i = yp;i??yp. Y p can be considered as a linear mapping from k-component
~? vectors to s-component y vectors. By regarding Y p as an observation operator of
the new k-dimensional space, we efiectively have linearized M. This linearization is
not the same as in the EKF equation (1.25). Rather, regarding Y p as an observa-
tion operator efiectively linearly interpolates between the full (possibly nonlinear)
observations yp;i.
Although the predicted ensemble may not be Gaussian, we apply the standard
Kalman fllter update equations (1.22) and (1.23) by making the substitutions
??p ! ?~?p = 0;
P ! ~P = (k?1)?1I;
H ! Y p;
H??p ! ?yp:
Application of these substitutions to equations (1.22) and (1.23) results in the EnKF
update equations for computing the updated covariance ~U and mean ?~?u.
~U = ?(k?1)I + (Y p)T R?1Y p??1; (1.37)
23
?~?u = ~U(Y p)T R?1(y ? ?yp): (1.38)
In principle, we could now compute the full state space mean and covariance
??u = ??p + Zp?~?u; (1.39)
U = Zp ~U(Zp)T: (1.40)
However, it would be an ine?cient to compute U directly; all that is required is to
obtain the updated ensemble ?u;i, and this can be obtained directly from ~U and ??u.
A Monte-Carlo method of obtaining the ensemble would be to randomly sam-
ple k points in the k-dimensional space according to a PDF with mean ?~?u and
covariance ~U. These points could them be transformed into the state space us-
ing (1.36), forming the updated ensemble. However, it is generally beneflcial if
the choice of updated ensemble is deterministic. We seek a linear transformation
from the predicted ensemble perturbations ??p;i (equivalently, the matrix Zp) to the
updated ensemble perturbations ??u;i (equivalently, the matrix Zu) such that the
covariance of the updated ensemble perturbations, as computed by (1.31), gives the
required U of equation (1.40). The transformation is written Zu = ZpW. Plugging
this into (1.31) gives U = (k? 1)?1Zu(Zu)T = (k? 1)?1ZpWW T (Zp)T . Thus if
W satisfles
~U = (k?1)?1WW T; (1.41)
equation (1.40) will be satisfled. There are several W matrices which satisfy equa-
tion (1.41). However there is a unique real W which satisfles (1.41) and is sym-
metric; this is our choice of W. It can be found numerically by computing the
24
Eigendecomposition of ~U,
~U = V DV T; (1.42)
leading to
W = pk?1V D1=2V T: (1.43)
The flnal updated ensemble is
?u;i = ??u +
kX
j=1
Wij??p;j; (1.44)
where the Wij are elements of the W matrix, and ??u is given by (1.39). This step
completes the EnKF update.
1.4.1 Parameter Estimation
There is a straightforward extension of the ensemble methods for cases in
which some model parameters are unknown. Consider the model
?j+1 = G(?j;p) (1.45)
and the extended state space to vectors having the form  =
2
66
4
?
p
3
77
5, where p, the
vector of model parameters, is treated as a state variable. The extended model
evolves as
 j+1 =
2
66
4
?j+1
pj+1
3
77
5 =
2
66
4
G(?j;pj)
pj
3
77
5 = ^G( j): (1.46)
Estimates of  (and therefore the parameters p) result from an implementation in
the same way as for ?, but in the space of  vectors.
25
1.4.2 Discussion
Although this form of EnKF has many advantages over the EKF, it may still
be computationally infeasible for systems exhibiting high dimensional chaos. As
the system size grows and the dynamical degrees of freedom increase, so too must
the number of ensemble members k in order to sample the state space. This is a
major drawback of ensemble methods like the EnKF, preventing their use for spa-
tiotemporal chaos in large domains which would require an unfeasibly large k. In
this case, the bottleneck in the computation may not be the EnKF update, but
rather the simulation of k difierent system states. This is especially true when the
model must simulate a system of partial difierential equations, as in many spa-
tiotemporally chaotic systems. This motivated the development of the LETKF of
Chapter 2. The LETKF algorithm is presented in the next section in a context
speciflc to Rayleigh-B?enard convection.
26
Chapter 2
Application of Data Assimilation to Rayleigh-B?enard Convection
The EnKF reduces the rank of covariance matrices by only keeping track of
the PDF in the space spanned by the ensemble members. The LETKF reduces the
rank further, by considering that, for spatiotemporal chaos with a small correlation
length, the system state at any point is uncorrelated with the system state at points
far away. For example, it is unnecessary to compute or store elements of covariance
matrices corresponding to the temperature at two points far from one another (we
can expect this element of the matrix to be approximately zero). Alternatively,
one can consider that in spatiotemporally chaotic systems the dynamics in local
regions tend to lie in a low dimensional space. To take advantage of this property
of spatiotemporal chaos, the LETKF performs local updates on overlapping patches
covering the domain. This is advantageous since the number of ensemble members
required is independent of the system size, making the method applicable to large
domains [10, 11]. Local regions may be spanned by fewer ensemble members, and
thus the required k is smaller for the LETKF as compared to the EnKF. Other than
this important difierence, the LETKF operates similarly to the EnKF. The LETKF
update uses the same equations as the EnKF, using an ensemble of system states to
represent the PDF. The only difierence being that the computations take place for
many local ensembles, rather than one global ensemble. The full LETKF algorithm
27
is presented in Section 2.1.
Our goal is to determine the full  uid state of a Rayleigh-B?enard convection
experiment, from a time series of shadowgraph measurements. We view this as a test
case investigation of the general usefulness of the LETKF technique for laboratory
experiments on spatiotemporal chaos. The  uid state ? consists of the variables  
and u deflned on the grid points (rm;`n;zl) of the cylindrical mesh; symbolically
? =
2
66
4
 
u
3
77
5:
For the model G we simulate the Boussinesq equations (1.1). We assume that mea-
surements are made at constant intervals ?t(tj ?j?t). In this context the measure-
ments come as a collection of shadowgraph pixels, y = [I(x1;y1) I(x2;y2) ::: I(xs;ys)]T
where I(xl;yl) is the light intensity at the location (xl;yl) of pixel l. Note that the lo-
cation of these intensity measurements need not occur on a uniform mesh; we assume
that their location is flxed but arbitrary. We assume for simplicity that R =  2I,
i.e., a multiple of the identity matrix, so that measurement noise is homogeneous
and uncorrelated with a standard deviation of  . This assumption is justifled for
most experimental setups, as shown in Section 3.2. The map M(?) outputs the
vector of pixel intensities y using a flnite resolution approximation to (1.6), where
for r2? we use a flnite difierence on the cylindrical mesh. Note that, since we require
kar2? k? 1 for (1.6) to be a good approximation, M is only weakly nonlinear.
An application of any of the previously introduced data assimilation techniques
(EKF, EnKF) to high aspect ratio Rayleigh-B?enard convection leads to enormous
computation. Either inversion of an N ?N matrix, or the requirement for a large
28
number of ensemble members make both methods infeasible. For example, we found
for the Rayleigh-B?enard problem (with ? = 20) that, using the EnKF, it was not
computationally feasible to use large enough ensembles to obtain results of any use.
For a small domain (? = 6 with ? = 2:0) we found that the EnKF required k ? 100
ensemble members. Since the dimensionality of Rayleigh-B?enard convection has
been shown to be extensive [20], we expect on the order of k ? 1000 ensemble
members would be required for an aspect ratio ? = 20.
2.1 The Local Ensemble Transform Kalman Filter
We now describe the LETKF?s update step (1.34). This Appendix is an adap-
tation to our Rayleigh-B?enard problem of the technique developed in [12]. Let ?mn
be a vector whose components consist of the collection of all elements of ? (in any
order) that lie on grid points within a horizontal distance L of the point (rm;`n) of
the mesh used by the model. We call ?mn a local state and L the local region radius.
There are as many (overlapping) local regions as horizontal grid points (rm;`n), see
Fig. 2.1. Note that, since the problem of interest is essentially two dimensional,
local regions are indexed by two indices (m;n). The three dimensional nature of the
system is re ected in the fact that, for each horizontal grid point, the vector ?mn
contains the state at all z levels zl; l = 1:::Nz. Associated with the updated and
predicted (global) ensemble members ?u;i and ?p;i are the local ensemble members
?u;imn and ?p;imn (all local states, global states, and ensemble states have an implied
time index j).
29
L
(r m , ?n )
Figure 2.1: Two local regions are shown on a reduced resolution mesh. Every grid
point (m;n) is the center of a local region. Associated with each local region (m;n)
is the local state vector ?mn consisting of state variables on the indicated horizontal
grid points and all vertical grid points associated with them.
Recal that the observation ensemble of predicted shadowgraphs fyp;1::: yp;kg
is deflned as yp;i = M(?p;i) (the projection of the predicted ensemble into the
observation space). Let yp;imn be all elements of yp;i within the local region (m;n).
If there are smn measurements made within the local region (m;n), then the vector
yp;imn has dimension smn. Form the matrix Y pmn ? [?yp;1mn j?yp;2mn j ??? j?yp;kmn] where
?yp;imn = yp;imn ? ?ypmn and ?ypmn is deflned as in (1.30). The local measurements ymn
have an associated local smn ?smn covariance matrix Rmn, which is equal to  2
multiplied by the smn ?smn identity matrix. We modify this matrix by forming the
30
tapered diagonal covariance matrix Qmn having i;ith element
Qiimn ? [ e(r=rf)2]2; (2.1)
where r is the (horizontal) distance from the grid point (m;n) to the location of
the ith element of ymn, and rf is some fallofi distance. This modiflcation efiec-
tively weights measurements further from the grid point (m;n) less heavily when
estimating the state of the (m;n) local region. This type of distance-dependant
modiflcation to the covariance matrices has been investigated previously [17]. The
form (2.1) in particular is very efiective for reducing the convergence time in our tests
on Rayleigh-B?enard convection when compared to other forms (for example a linear
function of r=rf). We also weight current measurements more heavily than prior
ones by the method of multiplicative variance in ation. The predicted covariance
matrix is in ated by a factor ?2 > 1, to lessen the in uence of prior measurements
on the current state, and to compensate in some rough way for model nonlineari-
ties [16, 12]. Ordinarily ? is chosen empirically.
We proceed to compute the updated ensemble. The following equations take
place in the k-dimensional space of local predicted perturbations ??p;imn. The inputs
are the global predicted ensemble ?p;i and the measurement y. The output is the
global updated ensemble ?u;i.
Compute each yp;i = M(?p;i).
For each grid point (m;n):
31
1. Form ymn from the elements of the measurement y, along with the tapered
covariance matrix Qmn.
2. Form ?ypmn and Y pmn from the ensemble yp;i.
3. Compute the updated k?k covariance matrix,
~Umn = [(k?1)??2I + (Y pmn)TQ?1mnY pmn]?1: (2.2)
4. Next compute the update vector
wmn = ~Umn(Y pmn)TQ?1mn(ymn ? ?ypmn); (2.3)
which is a transformation of the difierence between the predicted and actual
measurement into the k dimensional space of perturbations ??p;im;n.
5. Calculate the matrix
W mn = [(k?1) ~Umn]1=2 + wmn; (2.4)
where, by adding a vector to a matrix we mean adding it to each column of
the matrix. The 1=2 power here indicates taking a symmetric square root.
6. Form the local predicted ensemble f?p;1mn::: ?p;kmng from the global predicted en-
semble f?p;1::: ?p;kg as well as the matrix Zpmn ? [??p;1mn j??p;2mn j ??? j??p;kmn].
7. Finally, compute the local updated ensemble perturbations ??u;imn,
Zumn = ZpmnW mn: (2.5)
32
As before, Zumn ? [??u;1mn j??u;2mn j ??? j??u;kmn], and the local updated ensemble
is given by ?u;imn = ??pmn +??u;imn.
To complete the update step, the global updated ensemble member ?u;i is
obtained from the center elements of each ?u;imn. by setting the value of ?u;i at each
point (m;n) equal to the elements of ?u;imn corresponding to the center of local region
(m;n) (i.e., corresponding to the local state at grid point (m;n)). Note that each
local region is assimilated independently, allowing for massive parallelization.
To estimate parameters, simply replace ? with  everywhere in the above
steps. This formulation assumes state variables are spatially extended. Thus,
when adding global parameters to the state space we must assume that they are
spatially dependant. That is, when estimating both the Rayleigh number and
a of (1.6) (p = [R a]T), the state  is the concatenation of ? and ^p, where
^p = [R11 ::: Rmn ::: a11 ::: amn :::]T. The LETKF is then augmented by av-
eraging these parameter values over the grid, after the update step, to form global
parameters. This average is performed for each global ensemble member  u;i by
setting its ^p component to ^pu;i = [ ?Ri ?Ri ::: ?ai ?ai :::]T, where ?Ri and ?ai are the
ith ensemble members spatial averages of R and a, respectively. If the model allows
for spatially dependent parameters then this last averaging step is not necessary.
33
2.2 Direct Insertion
In order to assess how well the LETKF method is performing, we will compare
it to a more naive approach that we call Direct Insertion (DI). The DI approach is
motivated by a style of synchronization used by Pecora and Carroll in [35] in which
the available measurements simply substitute their corresponding state variables.
This only applies to when state variables are measured directly. With shadowgraph
measurements, no state variables are measured directly; however, there is a one to
one correspondence between a shadowgraph and the vertically averaged fleld ? (x;y).
With this in mind, the DI update step adjusts the predicted t = tj temperature fleld
? pj (x;y) to re ect the current measurement exactly.
At the time tj of the shadowgraph measurement Ij(x;y), the DI method up-
dates the predicted temperature fleld  pj (x;y;z) by adding a correction ? j(x;y;z)
which is the unique fleld that is quadratic in z, matches the boundary conditions at
jzj = 12, and for which the updated fleld  uj (x;y;z) =  pj (x;y;z)+? j(x;y;z) satisfles
Ij(x;y) = I?(x;y)1?ar2
?? 
u
j (x;y)
: (2.6)
This gives the update
? j(x;y;z) = (? uj (x;y)? ? pj (x;y))
 3
2 ?6z
2
?
; (2.7)
where ? uj (x;y) is found by solving a Poisson equation,
r2? uj (x;y) = 1a
?
1? I?(xc;yc)I
j(xc;yc)
?
; (2.8)
and (xc;yc) is the location of the closest pixel to (x;y) that is observed. Note that
with DI the velocity is not updated, uuj (x;y;z) = upj(x;y;z), rather it develops
34
through coupling with the temperature during the simulation step,
2
66
4
 pj+1(x;y;z)
upj+1(x;y;z)
3
77
5 = G
0
BB
@
2
66
4
 uj (x;y;z)
uuj (x;y;z)
3
77
5
1
CC
A:
The z-dependence of the predicted fleld is only slightly afiected by the update
since, if measurements are su?ciently frequent, the correction ? j(x;y;z) (which
has quadratic dependence) is small. This method is the most successful data assim-
ilation technique we have tested that does not use an update based on the Kalman
Filter. It is meant to represent what one might try without knowledge of the local-
ization techniques described in this section.
35
Chapter 3
Results
3.1 Perfect Model
In this section we describe so-called perfect model tests in which a time series
of states, generated from a Boussinesq simulation (? = 20, ? = 1, Pr = 1) of one
particular initial condition, serves as a proxy for the experimental system. Simulated
shadowgraph measurements of this time series are generated every ?t = 1=4 by
using (1.6) with the parameters a = 0:08, I?(x;y) = 0:5. By this technique we
generate a situation in which the ?true state? to be estimated and the model used to
estimate it both evolve under exactly the same dynamical rules. In Sec. 3.2 we use
real (not simulated) observations of a physical system for which the model dynamics
is surely not an exact description.
To reproduce the efiects of measurement noise we add to each pixel a small
random error that is an uncorrelated normally distributed number with mean zero
and standard deviation  . Measurements are made sparse by removing shadow-
graph pixels, leaving only those which lie on observation locations. We introduce
the measurement density ? ? s=(??2) as a measure of sparseness, where s is the
number of observation locations. For ? ? 4 we randomly and uniformly distribute
observation locations over the disk, otherwise the observation locations are placed on
a Cartesian grid covering the disk (giving more repeatable results when using sparse
36
measurements). The observation locations are flxed for the entire data assimilation
process.
We apply the LETKF and DI methods to our simulated shadowgraph time se-
ries to approximately reconstruct the original time series of true states. Performance
is quantifled via the temperature and mean  ow RMS relative error,
E (t) =
s
hj (x;y;z;t)? t(x;y;z;t)j2i
hj t(x;y;z;t)j2i ; (3.1)
E?u(t) =
s
hj?u(x;y;t)? ?ut(x;y;t)j2i
hj?ut(x;y;t)j2i ; (3.2)
where  t(x;y;z;t) and ?ut(x;y;t) are the temperature and mean  ow from the ?true?
time series of states, and h?i indicates a spatial average.
Simulated shadowgraphs are assimilated at times tj, j = 1:::J. During this
process we converge on an estimate of the system state (J chosen large enough to
ensure convergence). At time tJ assimilation is turned ofi and the flnal updated
state estimate is used as an initial condition for a long term forecast. We note
that the initial state estimate does not attain the minimum error, instead it occurs
about 1 tv into the forecast. This is a result of the simulation rapidly balancing the
flelds by strongly suppressing fleld errors with large wave numbers. This efiect is
very slight in the LETKF forecasts, but can be quite strong in DI forecasts. Three
measures of the quality of a state estimate are used: the minimum values attained
by E (t) and E?u(t) during a forecast, denoted Emin and Emin?u , and the predictability
time ?, deflned as the time when E (t) flrst crosses the (somewhat arbitrary) value
of 0:15. The perfect model tests reported here for the LETKF use k = 18 ensemble
members, a variance in ation factor ? = 1:0 ? 1:1, a local region radius L = 2:6d,
37
and a fallofi distance rf = 1:4d.
3.1.1 Performance with Noise/Sparseness
Here we document the performance of DI and the LETKF as a function of
measurement noise  and measurement density ?. We deflne the ideal scenario as
measuring a shadowgraph every tv=4 with ? = 127 (corresponding to a 451 ? 451
shadowgraph image) and  = 0:01 (this situation can be achieved in an experiment).
Under these conditions the DI and LETKF (with k = 18 ensemble members) con-
verge on a state estimate within ? tv and ? 3tv, respectively (observing ? 4 and
? 12 shadowgraphs, respectively). Both DI and the LETKF are efiective for es-
timation of the (unobserved) mean  ow ?u(x;y), however the LETKF achieves a
minimum error Emin?u that is less than half that of DI. The forecast error for a typical
state estimate is shown in Figs. 3.1 and 3.2. The general character of the forecasts
is a shadowing of the true state, followed by rapid divergence. When divergence
begins, the spatial structure of the error is concentrated near defects. This behavior
is expected, as the magnitude of the Lyapunov vector associated with the largest
Lyapunov exponent is largest at the location of defects [20].
38
 0
 0.2
 0.4
 0.6
 0.8
 1
 0  10  20  30  40  50  60  70  80
E
?
(t) forecast
t (t v )
DI
LETKF
 0  1  2  3
 0
 0.1
 0.2
t (t v )
DI
LETKF
E ?(t) converging
Figure 3.1: The error of the forecast temperature E (t) with  = 0:01 and ? =
127. Also shown in the small graph is E (t) as each method converges on a state.
Assimilation is turned ofi at time tJ = 3:25 in the small graph, corresponding to
time 0 in the large graph. The dashed line is our chosen threshold, E (t) 6 0:15,
below which we consider the forecasts ?good?.
39
 0
 0.2
 0.4
 0.6
 0.8
 1
 0  10  20  30  40  50  60  70  80
E
u
-
(t) forecast
t (t v )
DI
LETKF
 0  1  2  3
 0
 0.1
 0.2
t (t v )
DI
LETKF
E u-(t) converging
Figure 3.2: Mean  ow error E?u(t) of forecasts with  = 0:01 and ? = 127. Assimi-
lation is turned ofi at time tJ = 3:25 in the small graph, corresponding to time 0 in
the large graph. The insert shows E?u(t) as each method converges on a state.
Under non-ideal conditions the LETKF proves much more robust than DI.
Results for sparse measurements, shown in Fig. 3.3, demonstrate the large range of
?for which the LETKF converges. One can observe the existence of a critical density
of observations above which the LETKF does not substantially improve and below
which it fails to converge. By adjusting the parameters of the LETKF?s update step
(as described in the Appendix) we have been able to push the critical density as
low as ? = 1:3 without a signiflcant loss of quality in the state estimate. DI on the
other hand exhibits a steady increase in Emin and Emin?u as ? is decreased, as well as
a rapidly deteriorating forecast when even a few observation locations are removed.
40
 0
 20
 40
 60
 0 20 40 60 80 100 120
?
?
DI
LETKF
0.00
0.03
0.06
0.09
E
min u DI
LETKF
0.00
0.01
0.02
0.03
0.04
0.05
E
min ? DI
LETKF
Figure 3.3: Emin , Emin?u , and the predictability time ? as the density of observations
? is reduced.
Just as there is a critical measurement density, we have also found evidence
of a critical measurement frequency. This frequency lies somewhere between 1 and
2 shadowgraph images per vertical difiusion time for repeatable convergence of the
LETKF under ideal conditions. This corresponds to about 1 Hz in a typical exper-
iment.
Since the shadowgraph signal is simply variations from the background in-
tensity I?(x;y), the magnitude of measurement noise is best represented, not when
compared with the typical shadowgraph magnitude, but when compared to the typ-
ical RMS intensity variation  sg of a shadowgraph. In other words, the meaningful
signal to noise ratio is  sg= ( sg ? 0:12 when a = 0:08 and I?(x;y) = 0:5). DI relies
on the Poisson solve (2.8) which is fundamentally insensitive to noise (it smoothes
41
the right hand side). However, this insensitivity competes with the sensitivity of the
chaotic system dynamics when producing forecasts. The net result, in Fig. 3.4 indi-
cates that DI forecasts are only useful for  . 0:4 sg, whereas the LETKF operates
up to and exceeding  =  sg.
 0
 20
 40
 60
 0.2  0.4  0.6  0.8  1  1.2  1.4  1.6
?
?/ ?sg
DI
LETKF
0.00
0.03
0.06
0.09
E
min u
DI
LETKF
0.00
0.01
0.02
0.03
0.04
E
min ?
DI
LETKF
Figure 3.4: Emin , Emin?u , and the predictability time ? as measurement noise is in-
creased.
We note that all results are from one particular realization of the possible
?true? time series, generated from one particular initial condition. These results are
typical of what one can expect; however, variability can be expected (particularly
in ?) for difierent data sets.
42
3.1.2 Parameter Estimation
The relevant parameters available for estimation include not only the model
parameters Pr and R but also the observation operator parameter a of (1.6). In
general, the LETKF facilitates estimation of observation operator parameters in
exactly the same way as model parameters, by replacing M(?) by M(?;p) ? ^M( ).
The initial ensemblef p;10 ::: p;k0 gis constructed as before, from states sampled from
the attractor in the ? component, while the p component is sampled from a normal
distribution. When estimating a and R (with true values a = 0:08 and R = 3414)
the initial distribution was given mean (a = 0:07;R = 3073) and standard deviation
( a = 0:02; R = 683). The convergence process is demonstrated in Fig. 3.5.
43
 3310
 3350
 3390
 3430
 3470
 0  1  2  3  4  5  6  7  8
R
t (t v )
R
 0.073
 0.075
 0.077
 0.079
 0.081
 0  1  2  3  4  5  6  7  8
a
t (t v )
a
Figure 3.5: Simultaneously estimating the parameters a (with true value 0.08) and
R (with true value 3414.0). The error bars give a visual representation of the
ensemble spread, extending one standard deviation up and down. The thick error
bars represent the case ? = 127 and  = 0:01, while the thinner represent the sparse
measurement case, ? = 3:6 and  = 0:01.
Under ideal conditions the ensemble converges in 8tv on p = [R a]T =
[3414:26 0:07979]T?[1:61 0:000072]T, compared to the true value p = [3414:0 0:08]T.
Here the error estimates for R and a are the standard deviations of the p component
of the ensemble after the last update. Remarkably, even when measurements are
44
sparse (? = 3:6, near the critical measurement density) the parameter estimates are
very good, p = [3416:71 0:07976]T ? [9:9 0:00044]T. When estimating the state
and parameters simultaneously, the values of Emin and Emin?u are similar to those
shown in Fig. 3.3 and 3.4. That is, the ability to estimate the system state is not
adversely afiected when parameters are simultaneously estimated. It is important
to note that estimating parameters (in  space) requires more ensemble members
than when parameters are known, thus parameter estimation tests were performed
with k = 20.
Due to limitations in our simulation, we were unable to estimate Pr in the
full ? = 20 system. Instead we estimated Pr, a, and R together in small aspect
ratio tests (? = 4, R = 8540, Pr = 1). This small aspect ratio made possible
the application of the EnKF, allowing for comparisons between the LETKF and
the full EKF. Fig. 3.6 and 3.7 show the convergence of the ensemble toward the
true parameter values. Both the EKF and LETKF assimilated data from simulated
shadowgraphs (? = 93,  = 0:01), and both achieve lower than 1% error in Pr
parameter estimates.
45
2.5 5 7.5 10 12.5 15 17.5
t
8400
8450
8500
8550
8600
8650
8700
R
EKF
LETKF
Figure 3.6: Convergence of the ensemble mean of the parameter R toward the true
value, indicated by the red line. The initial ensemble had a large spread about the
mean R = 7857.
2.5 5 7.5 10 12.5 15 17.5
t
0.99
1
1.01
1.02
R
EKF
LETKF
Figure 3.7: Convergence of the ensemble mean of the parameter Pr toward the true
value, indicated by the red line. The initial ensemble had a large spread about the
mean Pr = 1:2.
46
3.2 Experiment
3.2.1 Experimental Setup
The experiment difiers from a perfect model scenario in that G and M are now
approximations, requiring robustness to model error as well as observation operator
error. In particular, the Boussinesq model is an approximation to the more exact
Navier-Stokes equations and our geometric optics treatment is an approximation to
a more involved physical optics treatment. For example, the Boussinesq equations
do not treat the temperature dependence of the  uid viscosity, thermal expansion
coe?cient, speciflc heat (at constant pressure), or conductivity; each of which varies
by 5% to 10% over the temperature range ?T of the experiment.
The geometry, parameter values and boundary conditions are closely matched
between experiments and simulations. For our experiments, the  uid is a thin
(d = 0:0602 cm) layer of carbon dioxide gas compressed at a gauge pressure 31.58
bar. The layer is surrounded by a circular boundary of radius 1.25 cm. In the exper-
iment, the top, bottom and lateral boundaries are composed of sapphire, aluminum
and polyethersulfone, respectively; the thermal conductivities of the boundaries ex-
ceed that of the gas by at least an order of magnitude. For this  uid, the critical
temperature difierence for convection onset is ?T = 6:02 ?C and tv = 1:66 s. A
flxed temperature difierence ?T = 10:23 ?0:09 ?C is imposed across the layer at a
flxed mean temperature of 22.6 ?0:1 ?C. These conditions correspond to R = 2902
(? = 0:7), Pr = 0:97, and ? = 20:8.
DI and the LETKF were used to assimilate shadowgraph images from the
47
experiment. Images were taken every ?t = tv=5 (3.0 Hz) as 395?395 bitmaps (? =
90) having  sg = 0:06, a signal to noise ratio of approximately 20 ( = sg ? 0:05).
In experiments, the true  uid state is not available for directly ascertaining the
accuracy of state estimates. Instead, we generate a forecast of the state estimate,
using the model, and compare the predicted shadowgraph sequence to subsequent
measurements (M is applied to the forecast state every ?t). Shadowgraphs are
flrst flltered by removing high frequency components (wavelengths less than d=2).
We then threshold the image such that half the pixels are set to 1 (the remaining
half are 0). This flltering/threshold procedure is applied to both the predicted and
measured shadowgraph time series. The natural error measure is then the fraction
of pixels incorrectly predicted, denoted EI.
The results reported here using experimental data use an in ation factor of
? ? 1:4. The large variance in ation is used to account for some of the model error,
and greatly improves stability. In this section, in which we assimilate experimental
data, the LETKF uses the parameter values k = 18, L = 2:6d, and rf = 1:4d.
3.2.2 Analysis of Measurement Noise
Recall that the s component vector ? represents measurement noise (yj =
M(?j) + ?). It is is a random variable which, in order to apply the Kalman fllter
methodology, is assumed to be normally distributed with mean 0. Two shadowgraph
images (I(1)? and I(2)? ) were taken below onset (?< 0) and compared to estimate the
distribution of noise in shadowgraph images. The the shadowgraph light intensity
48
at the location of pixel [i;j] in each image is denoted I(1)? [i;j] and I(2)? [i;j], and
the average is computed as Iavg? = (I(1)? +I(2)? )=2. The distribution of the quantity
?ij = I(1)? [i;j] ?Iavg? [i;j] over the entire image is denoted f(?ij). Fig. 3.8 shows
the distribution f(?ij) is symmetric and normally distributed with mean zero and
standard deviation  ? 0:003. Thus the assumption made regarding the normality
of measurement noise is a good approximation.
-0.01 -0.005 0 0.005 0.01
dij
0
0.05
0.1
0.15
0.2
0.25
f
H
d
ij
L
Figure 3.8: Probability distribution of shadowgraph pixel noise ?ij.
The noise at each pixel is weakly correlated to noise at nearby pixels. This
correlation can be measured via the quantity
C(m;n) ?
P
ij
?
I(1)? [i+m;j +n]?Iavg? [i+m;j +n]
??
I(1)? [i;j]?Iavg? [i;j]
?
P
ij
?
I(1)? [i;j]?Iavg? [i;j]
?2 :
The one dimensional correlation function is denoted C(n),
C(n) ? 14
X
i2+j2=n2
C(i;j): (3.3)
Fig. 3.9 shows this correlation versus pixel distance.
49
0 1 2 3
n
0
0.25
0.5
0.75
1
C
H
n
L
Figure 3.9: Correlation of shadowgraph pixel noise to the noise of nearby pixels.
There is a slight correlation between immediately adjacent pixels. Although we
could build this slight correlation into the noise covariance matrix R, Fig. 3.9 shows
that the correlation is small. Thus we assume that pixel noise is independent and
isotropic so that R is a multiple of the identity matrix. The full correlation C(m;n)
is plotted in Fig. 3.10. One can notice a slight asymmetry, the vertically adjacent
pixels are more highly correlated than horizontally adjacent ones. Presumably this
is speciflc to our shadowgraph apparatus.
50
-3
-2
-1
0
1
2
3
m
-3
-2
-1
0
1
2
3
n
100
10-1
10-2
10-3
CHm,nL
Figure 3.10: The correlation function C(m;n) of shadowgraph measurement noise
from the experiment.
The optimal value of  for the LETKF update step was typically much larger
than the measured value  = 0:003. The best results were obtained with  ?
0:01 ? 0:03, most likely due to the ability of an in ated  to roughly compensate
for some observation operator error.
3.2.3 Parameter Estimation and Performance with Sparseness
The LETKF is given 4tv to converge on state and R estimates, this is su?-
cient for both ideal (? = 90) and sparse observation (? = 4) cases. In the ideal
case, the LETKF converges on the parameter estimate R = 2625 (the experimen-
51
tally measured value is R = 2902 ? 26). When ? = 4 the LETKF converges on
the estimate R = 2491. These estimates are obtained consistently throughout the
experimental data set. Typical forecasts demonstrate a linear forecast error growth
up to the saturation point near EI = 0:5. Since R has been measured in the ex-
periment, parameter estimation is not necessary. However, forcing the LETKF to
use the measured R value harms the forecast, bringing it up to the level of the DI
forecast. This indicates that the advantage of the LETKF in this case lies in its
ability to estimate parameter values, as model error can typically be compensated
for, to some extent, by adjustment of model parameters ofi their measured values.
An example of the resulting mean  ow estimate is shown in Fig. 3.11. This
data assimilation technique has allowed us to obtain (indirectly) the mean  ow
from shadowgraph measurements. Fig. 3.12 shows a typical state estimate from
the LETKF. To the eye, DI state estimates look nearly identical to the LETKF
estimates. However, DI forecasts shown in Fig. 3.13, which use the ?true? value of
R = 2902 are signiflcantly worse than LETKF forecasts, which use their respective
R estimates.
52
Figure 3.11: Close up image of the mean  ow structure from the LETKF state
estimate inferred from experimental data. The background shows the experimental
shadowgraph image.
53
Figure 3.12: An estimate of the  uid state after assimilating for 4tv of data (J = 20
frames). A: The t = tJ shadowgraph measurement indicating columns of warm
rising  uid (dark) and cold descending  uid (light). B: Temperature proflle ? (x;y)
from the state estimate. C: The modeled shadowgraph M[ (x;y;z)] of the state
estimate for comparison to A. D: The inferred vorticity potential `(x;y) which
solves r2`(x;y) = ?^z ? (r? ?u) and indicates regions of clockwise rotating (dark)
and counter-clockwise rotating (light) mean  ow.
54
 0
 0.1
 0.2
 0.3
 0.4
 0.5
 0  5  10  15  20  25  30
E
I
t
? = 90
? = 4
DI
LETKF
Figure 3.13: Forecast error EI for DI and LETKF methods are shown for high
and low measurement densities. The LETKF forecast uses its parameter estimate
(R = 2625 for ? = 90,R = 2491 for ? = 4) while the DI forecast uses the measured
value R = 2902.
Fig. 3.14 shows how the forecasts of Fig. 3.13 compare with typical perfect
model forecasts (the best that can be expected) using the same parameters as the
experiment (R = 2902, ? = 20:8, Pr = 0:97) as well as the same measurement
frequency (?t = tv=5), density (? = 90), and noise level ( = sg ? 0:05). The
experiment?s forecast is unable to shadow the true state as seen in the perfect model
case, likely due to model error.
55
 0
 0.1
 0.2
 0.3
 0.4
 0.5
 0  10  20  30  40  50  60  70  80
E
I
t
DI (E) LETKF (E)
DI (PM)
LETKF (PM)
Figure 3.14: Forecast error EI for DI and LETKF methods in perfect model (PM)
tests and when using experimental data (E). The LETKF forecasts use the estimated
value of R, while DI uses the true value.
3.3 Discussion
The results from applying the LETKF and DI methods to experimental data
are far from the ideal results presented in the perfect model tests. There is signiflcant
model error. We have attempted to compensate for this model error by introducing
a temperature dependence in the thermal conductivity and speciflc heat (each varies
by about 5% over the temperature range in the experiment). However, we found
that this did not improve forecasts. The model error may be a result of other
assumptions made in the Boussinesq model. For example, that the velocity fleld is
divergence free.
Recently it has been shown [36] that the dimension density ?D = D=?2 (where
D is the dimension of the full system) for spiral defect chaos is D=?2 ? 0:25 when
56
? = 2:5. This result was obtained in the same disk-shaped geometry we have inves-
tigated, with aspect ratios up to ? = 15. This indicates that for our aspect ratio
of ? = 20, there is approximately D = 100 degrees of freedom (note that this is for
? = 2:5 only; the dependence of dimension on ? remains unknown). The number of
ensemble members will presumably scale with the number of dynamical degrees of
freedom in a local region, which can be computed from this dimension density as
Dlocal = ?DL2. We used a local region radius of L = 2:6d and a fallofi distance of
rf = 1:4d for the perfect model tests. For the experimental data, we found that a
local region radius of L = 2:6d and a fallofi distance of rf = 1:0d worked well. In
both cases, L is comparable to the correlation length of spiral defect chaos of 2:7d
when ? = 0:7 and 2:3d when ? = 1:0 [37]. With these local region radii, local regions
have (on average) only a few (1 or 2) degrees of freedom.
As the ensemble converges, it tends to conflne itself to a space of dimension
lower than k. The E dimension [38] DE is a measure for the number of important
directions in the space spanned by the ensemble perturbations. Roughly speaking,
it gives the dimension of the space in which the ensemble is spread out. The E
dimension is computed by forming the N?k matrix Z having the (global) ensemble
perturbations as its columns. The real eigenvalues of the k?k positive semi-deflnite
matrix ZT Z, denoted  2i are used to compute the E dimension [38],
DE ?
? kX
i=1
 i
!2? kX
i=1
 2i
!?1
: (3.4)
Fig. 3.16 and 3.15 show the E dimension decreasing as the LETKF converges on
the state and parameters. The update step generally increases the E dimension,
57
while the predict step generally decreases it. In perfect model tests the E dimension
decreases to approximately DE ? 5, while when working with experimental data, it
generally decreases to DE ? 8. This is consistent with the claim that local regions
contain only a few degrees of freedom, and indicates that one could optimize by
?pruning? the ensemble size as it converges. All the results in sections 3.1 and 3.2
are for a constantk = 18 (ork = 20 when estimating parameters), but we have found
that starting with k = 18 and reducing to k = 8 linearly within 10 measurement
times gives similar results with a signiflcant reduction in computation time. A large
number of ensemble members (k ? 20, as in the reported results of the previous
sections) are required only for the flrst few assimilation steps. In addition, the
strength of the model nonlinearities is largest when the ensemble spread is large
(during the flrst few assimilation steps) thus one might begin assimilation with a
large ? and reduce it linearly to speed convergence. This procedure was found to
be successful, but was not performed in the results reported here.
58
0 1 2 3 4 5 6
t
5
7.5
10
12.5
15
17.5
20
D
E
Figure 3.15: E dimension decreasing as the LETKF converges on perfect model
data. Blue dots indicate the time just after an update step, while the red dots
indicate the time just after the predict step. The ensemble has k = 20 ensemble
members, and the measurements are ideal (in the sense described in section 3.1).
59
0 0.5 1 1.5 2 2.5
t
5
7.5
10
12.5
15
17.5
20
D
E
Figure 3.16: E dimension decreasing as the LETKF converges on experimental data.
Blue dots indicate the time just after an update step, while the red dots indicate
the time just after the predict step. The ensemble has k = 22 ensemble members
in this example, and the measurements are at the highest measurement density (in
the sense described in section 3.2).
In our model for incomplete measurements, observation locations are dis-
tributed uniformly. However, data may be incomplete due to large data voids in
the shadowgraph time series. We have found that the LETKF can estimate the
 uid state in these regions when the data voids have a characteristic radius not
much larger than a correlation length. Fig. 3.17 shows a sequence of images of
the midplane temperature estimate taken from the LETKF?s assimilation of several
experimental shadowgraph images having a large data void in the upper right.
60
Figure 3.17: The midplane temperature fleld as the LETKF converges on an es-
timate of the system state using measurements which have a large data void. No
measurements occur in the region with a red tint. The sequence of images occur at
the times tj = t1; t3; t10; and t20.
We have investigated two methods for estimating the  uid state in Rayleigh-
B?enard convection experiments, DI and the LETKF. Both methods are efiective
for this purpose, with the LETKF outperforming DI both when using experimental
data and in perfect model tests, especially when data is sparse/noisy. The LETKF
is a promising technique for large experimental systems, as its complexity does not
grow with the system size. In addition, it can take advantage of multiple processors,
even if the model cannot, by parallelizing over the ensemble members during the
61
forecast step and over grid points during the update step. The LETKF method we
have presented is potentially applicable to a large class of spatiotemporally chaotic
laboratory experiments. Support for part I of this thesis we provided by National
Science Foundation Grant 04-34193 and 04-34193 and the O?ce of Naval Research
(Physics).
62
Part II
THE EFFECT OF ROTATION ON DYNAMO ACTION IN
MAGNETOHYDRODYNAMIC TURBULENCE
63
Chapter 4
Introduction to Magnetohydrodynamics
Magnetohydrodynamics (MHD) is concerned with the study of conducting
 uids (liquid metals or plasma), and the interaction between the  uid  ow and
electromagnetic flelds. One motivation for the study of MHD is that planets and
stars contain complex  ows of molten metal or plasma. These  ows are thought to
be driven by convection, and may be responsible for the large scale magnetization
of the atmosphere of stars and some planets such as the Earth. However, though
the fundamental physics of MHD is understood, there is still much to learn about
the dynamics of the  ows, which are often chaotic and complex. For example, much
research is focused on the question: Does an initial small magnetic fleld grow as a
result of the  ow? The answer to this question depends on the system investigated;
if the answer is yes, the system is known as a dynamo.
The study of MHD via computer simulations is generally separated into two
classes. In one case, the problem is simplifled as much as possible to bring out the
essential dynamics [39, 40, 41, 42, 43]. In these studies, to isolate the bulk dynamics
from the efiects of boundary conditions, periodic domains (x;y;z 2 [0;2?]) are
employed. Rotation and convection are usually not included in favor of external
forcing to drive the system. These studies investigate the role of turbulence and the
conditions necessary for the generation of a dynamo in the simplest situations. The
64
other case involves complex simulations in which realistic boundary conditions are
employed [44, 45, 46, 47]. These simulations generally include convection to drive
the system and, in some cases, attempt to model the entire Earth. The study of
MHD simulations (in both cases) is limited by the available resolution of computer
models, which in turn is limited by available computing power.
There is a region between the two cases of complex, fully modeled  ow and
simple, idealized  ow, which has not received as much attention as either extreme.
It is the aim of this chapter to investigate one aspect of the MHD problem in
this region which has not been addressed previously: the inclusion of the efiects of
rotation in an externally forced periodic domain. Investigation of this particular
problem is interesting, as the efiects of rotation may be isolated, providing insight
into otherwise complex phenomena. There is some evidence that rotation plays a
role in the ability of planets to generate a magnetic fleld. Venus, with a rotation rate
116 times slower than the Earth, has no measurable terrestrial magnetic fleld [48].
Mercury, with a rotation rate 176 times slower than earth, has a magnetic dipole
moment about 1700 times weaker than Earth?s [48]. We will focus on incompressible
motion in which the  uid mass density is constant.
The dimensionless incompressible MHD equations describing the evolution of
the velocity fleld v and magnetic fleld B take the form [49, 50]
@v
@t + v?rv = ?rp+ (r?B)?B +?r
2v + F; (4.1)
@B
@t + v?rB = B?rv +?r
2B; (4.2)
r?v = 0; (4.3)
65
r?B = 0: (4.4)
A constant external force F is present to drive the  ow,
F = [(sinz + cosy) ^x + (sinx+ cosz) ^y + (siny + cosx) ^z]: (4.5)
All quantities above are dimensionless. Our simulations are done on a peri-
odic domain in (x;y;z) where the periodicity length is 2?. Using overbars to denote
dimensional quantities, we have that (x;y;z) = (2?=?L)(?x;?y;?z) where ?Lis the dimen-
sional periodicity length of the system. Our dimensionless system with periodicity
length 2? lends itself to the spectral decomposition of the flelds described in Chap-
ter 5 in which the allowed wave numbers in each coordinate direction are integer
kx;y;z 2 N. The dimensionless force as deflned in equation (4.5) has a spatial RMS
value of p3. This implies that velocities are normalized to ?U ?
q
?L?F=(2?p3??),
(i.e., v = ?v=?U) where ?F is the spatial RMS of the dimensional applied force density
and ?? is the mass density of the  uid (assumed constant in space and time). Times
are in units of ?L=(2??U) (i.e. t = 2??U?t=?L) and the magnetic fleld is normalized to
p??? ?U (i.e., B = ?B=( ?Up???)), where ? is the magnetic permeability of the  uid.
In these units, a normalized magnetic fleld B with magnitude jBj = 1 represents a
magnetic fleld for which the dimensional Alfv?en wave velocity is ?U; thus these units
are termed Alfv?enic units.
The flrst three equations (4.1), (4.2), and (4.3) are all that are necessary to
simulate the  ow, since if the initial B fleld is divergence free, equation (4.2) implies
that B remains divergence free for all time. The dimensionless pressure p is found by
taking the divergence of equation (4.1), imposing (4.3), and solving the Poisson-type
66
equation
r2p = r?F?r?(v?rv) +r?((r?B)?B): (4.6)
The quantities ? = 2???=?U?L and ? = 2???=?U?L which appear in (4.1) and (4.2)
are the dimensionless viscosity and magnetic difiusivity respectively; ?? and ?? are the
dimensional viscosity and magnetic difiusivity respectively. The Reynolds number
and magnetic Reynolds number are
R =
?L?Urms
?? =
2?Urms
? and (4.7)
Rm =
?L?Urms
?? =
2?Urms
? : (4.8)
Here, ?Urms and Urms are the space/time RMS average of the velocity fleld ?v(?x;?t) and
v(x;t) respectively, Urms = ?Urms=?U. Larger R (smaller viscosity) leads to smaller
structure in the v  ow and generally larger  ow velocities. Similarly, larger Rm
(lower Ohmic resistance) leads to smaller structure in the B fleld and less dissipation.
The ratio Rm=R = ?=? is known as the magnetic Prandtl number Prm.
Dimensional analysis gives the viscous dissipative scale ??v = ??3=4=??1=4 and
Joule difiusive scale ??B = ??3=4=??1=4, where ?? is the input power per unit mass [39].
These scales are the characteristic size of the smallest structures in the ?v and ?B fleld,
respectively. The ratio of the viscous dissipative scale to the Joule difiusive scale
is ??v=??B = (??=??)3=4 = (?=?)3=4 = Pr3=4m . For liquid metals the magnetic Prandtl
number is very small [44]Prm ? 10?5?10?6, and thus structures in the  uid  ow are
much smaller than structures in the magnetic fleld. Such a large separation in scales
means that a fully realistic simulation must employ a very dense mesh, capable of
resolving the  uid  ow while simultaneously being large enough to capture the large
67
spatial scales of the magnetic structure.
The forcing (4.5) is an oft studied form [43, 42, 51]. In the absence of a
magnetic fleld, and if su?ciently weak, it generates a so-called ABC  ow (vABC =
F=?). This  ow has as its only wave numbers jkxj = 1, jkyj = 1, and jkzj = 1 and
is known to generate chaotic particle trajectories. The ABC  ow also possesses the
special property that r? vABC = 0, r? vABC = vABC, and thus (r? vABC) ?
vABC = 0. This  ow has the maximum possible helicity for flxed kinetic energy,
where helicity is deflned as
H ? v?(r?v): (4.9)
Flow helicity is thought by some to be important for dynamo generation of magnetic
flelds [52]. For this ABC  ow the helicity and energy density coincide. At su?ciently
small ?, the ABC  ow becomes unstable, and turbulent  ows with time-dependence
and higher wave number components result.
The dimensionless kinetic energy Ev and magnetic energy EB in the periodic
domain ? are
Ev =
Z
?
1
2jvj
2dx3; (4.10)
EB =
Z
?
1
2jBj
2dx3: (4.11)
Energy is injected by the external forcing F into the  ow v, after which it cascades
to smaller scales and dissipates through viscous damping or is transferred to the B
fleld. The energy in B in turn is damped through Ohmic dissipation.
68
4.1 The Kinematic Dynamo
We note that the system of equations (4.1)-(4.4) has a solution with the mag-
netic fleld identically equal to zero. In such a case the  uid velocity can (depending
on ?) come to a turbulent, chaotic, periodic, or time-independent steady state. One
can then ask whether this  ow state is stable to the introduction of a small initial
magnetic perturbation. If it is unstable, the initial perturbation will grow with time
and is expected to eventually saturate in a nonlinear magnetized state. The prob-
lem of determining the stability of an initially unmagnetized  ow to small magnetic
perturbations is referred to as the kinematic dynamo problem.
The magnetic difiusivity is critical in determining whether a  ow is a kinematic
dynamo. For ? << 1 the magnetic fleld is largely frozen into the  uid, stretching
through difierential  ow in v, while difiusion is unimportant in all but the smallest
scales (this situation is conducive to dynamo action). For ? >> 1 the magnetic
fleld is highly dissipative, decaying toward B = 0. Taking a dynamical systems
perspective, for ? >> 1 the system attractor lies in the B = 0 plane. A blowout
bifurcation occurs for a critical value of ? at which the attractor expands into the full
v; B space [42]. In order to study the growth of a small initial B fleld, we B-linearize
equations (4.1) and (4.2) around B = 0. Since the Lorentz term in equation (4.1)
is second order in B, we arrive at the standard Navier-Stokes equation,
@v
@t + v?rv = ?rp+?r
2v + F: (4.12)
In our B-linearized situation, this equation is completely decoupled from the induc-
tion equation (4.2), which remains unchanged since every term is linear in B. The
69
(simplifled) kinematic MHD equations consist of (4.12), (4.2), and (4.3).
Because the velocity fleld now evolves according to equation (4.12), it takes
all the usual properties common to standard unmagnetized  uids. In particular, it
becomes turbulent for su?ciently small ?; if the  uid is not rotating, the energy
spectrum is isotropic and follows the well known Kolmogorov spectrum Ev ?k?5=3
for turbulence in the inertial range (the range where inertial forces dominate both
viscous forces and the external force F). Because v is time-dependant we cannot
expect to flnd exactly time-exponential B solutions. However, we can still investigate
the average growth or decay of EB.
4.2 Rotation
The efiect of rotation on dynamo action was flrst investigated by Mofiatt
in [53, 54]. Here we investigate the efiect of rotation in the kinematic case, in which
the study of the velocity fleld reduces to the study of rotating  uids without regard
to the magnetic fleld. This is advantageous because much literature already exists
on rotating  uid turbulence [55, 56].
In a frame rotating with constant angular velocity ? = ?^z, the addition of
the Coriolis and centrifugal forces modifles equation (4.12),
@v
@t + v?rv = ?rP ?2??v +?r
2v + F: (4.13)
The rotation rate ? is dimensionless and takes Alfv?enic units, ? = ???L=(2??U).
Since the centrifugal force ???(??r) can be written as a gradient of a scalar, it
is incorporated into the new ?pressure? P ?p? (1=2)j??rj2. The addition of the
70
Coriolis force leads to two natural dimensionless numbers, the Rossby and Ekman
numbers,
Ro ?
?Urms
2?L?? =
Urms
4??; (4.14)
Ek ? ??2?L2 ?? = ?8?2? = RoR : (4.15)
The Rossby number measures the relative magnitude of inertial forces to the Coriolis
force, while the Ekman number measures the relative magnitude of viscous forces
to the Coriolis force.
When Ro ? 1 and Ek ? 1 the  uid may be regarded as rapidly rotating.
Rapidly rotating  ows do not follow the Kolmogorov energy spectrum. Energy
typically falls faster with wave number, Ev ? k?2 [56]. In addition, the energy
wave number spectrum is highly anisotropic, with much smaller energy in the kz
direction (along the axis of rotation). The Taylor-Proudman theorem states that,
for the most rapidly rotating  ows, this efiect can be so strong that the  uid motion
becomes essentially two dimensional, with derivatives in z approximately zero. To
see how this happens, we note that Ro ? 1 implies that, to a flrst approximation,
the inertial terms in (4.13) are small compared to the Coriolis force, while Ek ? 1
implies that the viscous force is similarly small. Thus, to a flrst approximation, we
ignore all terms except the Coriolis force and the pressure, which in steady state
must balance each other, 2? ? v = ?rP. Removing the pressure by taking the
curl gives (??r)v = ?@v@z = 0. Thus, v is approximately independent of z for large
?. Two dimensional  ows are known to exhibit an inverse energy cascade to larger
scales than the scale at which energy is introduced [57].
71
Linearizing equation (4.13) about the state v = 0 and taking a curl we obtain
@
@t(r?v) = 2??rv +?r?(r
2v): (4.16)
This equation supports traveling wave solutions known as inertial waves [54, 55].
Letting v = v?ei(k?r?!t) gives
(k ?v?)! = 2i(??k)v? ?i?k2(k ?v?); (4.17)
which leads to the dispersion relation ! = 2?(kz=k)?i?k2 and the condition v??k =
0. This relation shows that inertial waves cannot have angular frequencies above
2?. The quality factor Q?Re(!)=(2Im(!)) of these waves is given by
Q = kz8?2k
? ?
2?2?
??1
; (4.18)
where ? = 2?=k. This expression shows that the wave with the highest Q factor
has kx = ky = 0, kz = 1. The factor ? ?2?2?? is similar to the Ekman number (and
exactly equal to it for the largest scale ? = 2? waves). Equation 4.18 holds only
in the bulk (far from boundaries). In bounded  ows the quality factor is known to
have a difierent dependence on Ek than derived here [58].
72
Chapter 5
Numerical Approach
5.1 Spectral Representation
The flelds v = vx^x + vy ^y + vz ^z and B = Bx^x +By ^y +Bz ^z are expanded in
Fourier series in which kx; ky; kz are integer,
vx;y;z(x) =
1X
kx;ky;kz=?1
v(k)x;y;zeik?x; (5.1)
Bx;y;z(x) =
1X
kx;ky;kz=?1
B(k)x;y;zeik?x; (5.2)
v(k) = v(k)x ^x + v(k)y ^y + v(k)z ^z; (5.3)
B(k) = B(k)x ^x +B(k)y ^y +B(k)z ^z: (5.4)
The flelds v and B are real, thus the conditions v(?k)x;y;z = (v(k)x;y;z)? and B(?k)x;y;z =
(B(k)x;y;z)? must hold, where ? is the complex conjugate operation. The condition
v(0)x;y;z = 0 is enforced, so that the  uid has no net momentum. We also require
B(0)x;y;z = 0 since the k = 0 mode does not feel the efiects of Ohmic dissipation
(r2B = 0 for this mode), and is thus unrealistic. For convenience we deflne the
nonlinear terms
a = v?rv; (5.5)
g = B?rv; (5.6)
h = v?rB; (5.7)
73
which have associated Fourier transforms a(k), g(k), and h(k), deflned similarly to
equations (5.1 - 5.4). Taking the Fourier transform of the rotating kinematic MHD
equations (4.13) and (4.2) we get the system of ordinary difierential equations:
d
dtv
(k) = ?ikP(k) ?a(k) ?2??v(k) ??k2v(k) + F(k); (5.8)
d
dtB
(k) = g(k) ?h(k) ??k2B(k); (5.9)
one set of equations for every discrete k vector. It is straight-forward to compute
the pressure in spectral space,
P(k) = k?2?ik ?a(k) ?2??(ik ?v(k))?: (5.10)
In this representation, the flelds v, B, a, g, h and P together represent the system
state. Although v and B alone are enough to determine all these quantities, it is
more convenient to think of all 6 flelds collectively as the system state. The method
used is pseudospectral, the nonlinear terms a, g, and h are computed in physical
space, while all derivatives are computed in spectral space.
The Fourier space representation of all flelds is truncated to include only wave
numbers satisfying jkj<K, reducing the problem to only the largest, most impor-
tant scales. Truncating scales smaller than the viscous dissipative scale K ? 2?=?v
has little efiect of the dynamics, as the removed modes have negligible energy con-
tent. The sums in equation (5.1) and (5.2) then only sum over modes satisfying
jkj<K. The physical space representation of flelds consists of values on aN?N?N
Cartesian grid.
74
5.2 Time Stepping
In this section we discretize the system in time, forming a map which takes
the system state at the current and previous times, t = tj;tj?1;tj?2:::, and outputs
the system state at the next time t = tj+1 = tj +?t. The state at the current time
t = tj consists of vj, Bj, aj, gj, hj, and Pj. We use a 3rd order semi-implicit method
which mixes the explicit Adams-Bashforth method with the implicit Adams-Moulton
method. For the system @y=@t = M(y) the 3rd order Adams-Moulton method is
yj+1 = yj +?t?(5=12)M(yj+1) + (8=12)M(yj)?(1=12)M(yj?1)?; (5.11)
while the 3rd order Adams-Bashforth method is given by
yj+1 = yj +?t?(23=12)M(yj)?(16=12)M(yj?1) + (5=12)M(yj?2)?: (5.12)
The Adams-Moulton method requires that M be inverted to solve for yj+1, hence
it is used for the linear terms, while the Adams-Bashforth method is used for the
nonlinear terms. For the system @y=@t = Ly + NL(y) with linear term Ly and
nonlinear term NL(y) the mixed method is
yj+1 = (1? 512?tL)?1
?
(1 + 812?tL)yj ? 112?tLyj?1 +
?t
h23
12NL(yj)?
16
12NL(yj?1) +
5
12NL(yj?2)
i?
: (5.13)
Although the induction equation (4.2) cannot produce divergence in the B
fleld when solved exactly, all time stepping schemes (including implicit methods like
the mixed method above) are unstable in this regard, and will generate a divergence.
This instability is easy to overcome, by adding a divergence cleaning step after each
75
time step. Divergence cleaning can be accomplished in spectral space by subtracting
the projection of B(k) onto k,
B(k) ! B(k) ? ^k(^k ?B(k)): (5.14)
The flnal time stepping method takes all 6 flelds and their Fourier transforms
for the current and last two time steps, and outputs all flelds and their Fourier
transforms for the next time step:
Input : vj;Bj;aj;gj;hj;Pj;v(k)j ;B(k)j ;a(k)j ;g(k)j ;h(k)j ;P(k)j ;
vj?1;Bj?1;aj?1;gj?1;hj?1;Pj?1;v(k)j?1;B(k)j?1;a(k)j?1;g(k)j?1;h(k)j?1;P(k)j?1;
vj?2;Bj?2;aj?2;gj?2;hj?2;Pj?2;v(k)j?2;B(k)j?2;a(k)j?2;g(k)j?2;h(k)j?2;P(k)j?2:
Output : vj+1;Bj+1;aj+1;gj+1;hj+1;Pj+1;v(k)j+1;B(k)j+1;a(k)j+1;g(k)j+1;h(k)j+1;P(k)j+1:
The following steps take place in the order they are written. Below, the notation
Ff(x) is used to indicate the Fourier transform of the function f(x), while F?1f(k)
indicates the inverse transform.
Let R(k)x = (1? 812?t?k2)(v(k)x )j +?t
h 1
12?k
2(v(k)
x )j?1 +
2 812?(v(k)y )j ?2 112?(v(k)y )j?1 +
23
12(?ikxP
(k)
j ?(a
(k)
x )j)?
16
12(?ikxP
(k)
j?1 ?(a
(k)
x )j?1) +
5
12(?ikxP
(k)
j?2 ?(a
(k)
x )j?2) +Fx
i
:
Let R(k)y = (1? 812?t?k2)(v(k)y )j +?t
h 1
12?k
2(v(k)
y )j?1 ?
2 812?(v(k)x )j + 2 112?(v(k)x )j?1 +
23
12(?ikyP
(k)
j ?(a
(k)
y )j)?
16
12(?ikyP
(k)
j?1 ?(a
(k)
y )j?1) +
76
5
12(?ikyP
(k)
j?2 ?(a
(k)
y )j?2) +Fy
i
:
Let R(k)z = (1? 812?t?k2)(v(k)z )j +?t
h 1
12?k
2(v(k)
z )j?1 +
23
12(?ikzP
(k)
j ?(a
(k)
z )j)?
16
12(?ikzP
(k)
j?1 ?(a
(k)
z )j?1) +
5
12(?ikzP
(k)
j?2 ?(a
(k)
z )j?2) +Fz
i
:
Compute (v(k)x )j+1 = (1 +
5
12?t?k
2)R(k)x + 2 5
12?t?R
(k)
y
(1 + 512?t?k2)2 + (2 512?t?)2 : (5.15)
Compute (v(k)y )j+1 = (1 +
5
12?t?k
2)R(k)y ?2 5
12?t?R
(k)
x
(1 + 512?t?k2)2 + (2 512?t?)2 : (5.16)
Compute (v(k)z )j+1 = (1 + 512?t?k2)?1R(k)z : (5.17)
Compute B(k)j+1 = (1 + 512?t?k2)?1
?
(1? 812?t?k2)B(k)j +
?t
h 1
12?k
2B(k)
j?1 +
23
12(g
(k)
j ?h
(k)
j )?
16
12(g
(k)
j?1 ?h
(k)
j?1) +
5
12(g
(k)
j?2 ?h
(k)
j?2)
i?
: (5.18)
Assign B(k)j+1 ! B(k)j+1 ? ^k(^k ?B(k)j+1) (5.19)
Compute vj+1 = F?1v(k)j+1 and Bj+1 = F?1B(k)j+1: (5.20)
Let J(k)mn = ikm(v(k)n )j+1 8m;n2fx;y;zg: (5.21)
Let L(k)mn = ikm(B(k)n )j+1 8m;n2fx;y;zg: (5.22)
Compute Jmn = F?1J(k)mn and Lmn = F?1L(k)mn 8m;n2fx;y;zg: (5.23)
Compute (an)j+1 =
X
m2fx;y;zg
(vm)j+1Lmn 8n2fx;y;zg: (5.24)
Compute (gn)j+1 =
X
m2fx;y;zg
(Bm)j+1Lmn 8n2fx;y;zg: (5.25)
Compute (hn)j+1 =
X
m2fx;y;zg
(vm)j+1Jmn 8n2fx;y;zg: (5.26)
Compute a(k)j+1 = Faj+1; g(k)j+1 = Fgj+1; and h(k)j+1 = Fhj+1: (5.27)
Compute P(k)j+1 = k?2
?
ik ?a(k)j+1 ?2??(ik ?v(k)j+1)
?
: (5.28)
77
Compute Pj+1 = F?1P(k)j+1: (5.29)
For small ?t this system approximates the continuous equations. There are a total
of 25 inverse Fourier transforms and 9 forward Fourier transforms, for a total of
34 transforms per time step. The transform in step (5.29) is not required, but is
included for completeness. Each Fourier transform is accomplished by a fast Fourier
transform (FFT) using a publicly available library [59] which takes advantage of the
fact that the flelds in physical space are real, speeding up transforms by a factor of
roughly 2 over a complex calculation.
5.3 Parameters and Resolution Requirements
Ideally we wish to resolve all scales up to ?v and ?B. However, in practice
this is usually not necessary, as the energy content of modes drops quickly in k,
even for scales larger than ?v or ?B. The required K must be large enough to
resolve the inertial range, but need not increase to the point where ?v is resolved.
When the efiects of rapid rotation are included, the resolution requirements drop
further. By the Taylor-Proudman theorem, rapidly rotating  ows do not require
high resolution in the ^z direction (the axis of rotation). However, although some
speedup could be gained by exploiting these properties, a uniform resolution (N3)
grid was always employed. The highest possible resolved wave number is given by
the Nyquist frequency K = N=2. However, in many studies involving the simulation
of partial difierential equations in a periodic cube, the relation K = N=3 is used to
avoid aliasing issues. It was found that both cutofi wave numbers resulted in nearly
78
identical behavior for the runs performed, thus the condition K = N=2 was used for
simplicity (and because most investigated resolutions were a power of 2).
We investigate the range ? ? 10?1 ?10?2, ? ? 100 ?10?1, and rotation rates
up to ? = 8. The inverse of ? appears in the Reynolds number. Thus, we typically
refer to ? by its inverse (and similarly with ?). In particular, the values ??1 = 6:3,
18, and 36 were used for most runs, for which the resolution was N = 16, 32, and
64, respectively. All three values of ? are small enough that without rotation the
ABC  ow vABC is unstable. The Reynolds numbers, Ekman numbers, and Rossby
numbers of these  ows will be computed from Urms and discussed in Chapter 6.
In the dimensionless form employed, a time ?t = 1 is on the order of one
large eddy turnover time, the characteristic time for the large scales in the  ow to
transit the domain. The most basic requirement for the time step is then ?t? 1. To
resolve the process of advection we must also have the standard Courant-Fredericks-
Levy (CFL) condition ?t < ?x=Urms = 2?=(NUrms). Another requirement is that
inertial waves with frequency 2? must be well resolved. The numerical method is
3rd order in time, and thus the eventual time step is chosen at a point well within
the region in which convergence follows a ?t?3 law. Typically this results in a much
smaller time step than the CFL condition requires. For example, most ??1 = 18
runs used ?t = 0:002, while most ??1 = 36 runs used ?t = 0:0005, about one order
of magnitude smaller than the CFL condition requires and giving at least 500 time
steps per period (? = 2?=?) for the fastest rotation rate (? = 8).
79
Chapter 6
Simulation Results
Simulations of the  uid  ow velocity v begin with the initial condition v(x;y;z;0) =
0 and evolve for many large eddy turnover times (typically to t = 400 ?T). When
calculating time averages, the evolution up to the time t = 100 ?t? is disregarded to
ensure that transient efiects are no longer present and that the dynamics re ect mo-
tion on the attractor. The investigated rotation rates are ? = 0; 0:5; 1; 2; 4; 6; 8.
The notation
h?i =
s
1
(2?)3
Z
?
j?j2dxdydz (6.1)
is used for the spatial RMS average of a scalar or vector quantity ?; and the notation
hh?ii = 1T ?t
?
Z T
t?
h?idt (6.2)
is used for averages over space and time.
The average RMS velocity Urms = hhvii is shown for ??1 = 36, ??1 = 18, and
??1 = 6:3 in Tables 6.1, 6.2, and 6.3, respectively. This average velocity is used to
compute the Reynolds number, Rossby number, and Ekman number. The largest
investigated Reynolds number is R? 4?103.
80
? 0 1 2 4 6 8
Urms 4.63 5.28 5.55 8.58 12.89 16.78
R 1048 1194 1256 1941 2916 3796
Ro (?10?1) 1 7.37 4.20 2.21 3.42 3.36
Ek (?10?4) 1 6.17 3.34 1.14 1.17 0.89
Table 6.1: Urms, Reynolds number R, Rossby number Ro, and Ekman number Ek
for ??1 = 36 and various rotation rates.
? 0 0.5 1 2 4 6 8
Urms 4.57 4.48 4.70 6.17 9.76 13.01 15.47
R 517 506 532 698 1104 1471 1750
Ro (?10?1) 1 7.13 3.74 2.45 1.94 1.72 1.53
Ek (?10?4) 1 14.1 7.03 3.52 1.76 1.17 0.880
Table 6.2: Urms, Reynolds number R, Rossby number Ro, and Ekman number Ek
for ??1 = 18 and various rotation rates.
81
? 0 0.5 1 2 4 6 8
Urms 5.53 5.12 5.45 6.35 8.10 8.90 8.90
R 219 203 216 251 321 352 352
Ro (?10?1) 1 8.15 4.34 2.53 1.61 1.18 0.89
Ek (?10?4) 1 40.2 20.1 10.0 5.03 3.35 2.51
Table 6.3: Urms, Reynolds number R, Rossby number Ro, and Ekman number Ek
for ??1 = 6:3 and various rotation rates.
One can see that the conditions for rapid rotation (Ek ? 1 , Ro ? 1) are
satisfled for these  ows. Equation (4.18), derived in section 4.2 shows that the
quality factor for inertial waves goes like Ek?1. Thus, for the investigated rotation
rates, especially those for which Ek ? 10?4, we can expect to see inertial waves
lasting for long times.
6.1 Flow Characterization
Each simulation stores a time series of the average helicity hHi, input power
hF?vi, inertia hv?rvi, Coriolis force h2??vi, viscous force h?r2vi, pressure force
hrPi, as well as the energy Ev, and EB. Fig. 6.1 shows the time series of four force
strengths from equation (4.13) for various rotation rates. Fig. 6.2 shows time series
of the average helicity, velocity, and external power input. All quantities exhibit a
chaotic time evolution.
82
Figure 6.1: Time series (from t = 100 to t = 400) of the average inertia hv ?rvi, pressure force hrPi, viscous force h?r
2
vi,
and Coriolis force h2??vi for various rotation rates (? = 0, 2, 4, 6, and 8) with ?
?1
= 18. For reference, hFi =
p
3.
83
Figure 6.2: Time series (from t = 100 to t = 400) of the average velocity hvi, helicity hHi, and power hF?vi for various rotation
rates (? = 0, 2, 4, 6, and 8) with ?
?1
= 18.
84
A convenient means of characterizing these plots is through the RMS  uctua-
tion
hh?ii =
? 1
T ?t?
Z T
t?
(h?i?hh?ii)2 dt
?1=2
: (6.3)
Figures 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, and 6.9 show these quantities as a function of ?
for the cases ??1 = 36, ??1 = 18, and ??1 = 6:3. In each case the error bars do not
indicate uncertainty in the mean quantity (this uncertainty is of the order of the line
thickness). Rather, they extend vertically to represent the time variability of each
quantity (as computed by (6.3)). We note that, when ??1 = 6:3, the variability
hh?ii in all  ow variables drops to near zero for ? > 4. For these fast rotation
rates, the ??1 = 6:3  ow is constant in time. Evidently the ABC  ow is stabilized
by the addition of the Coriolis force. However, the  ow is turbulent and chaotic for
all investigated ? when ??1 = 18, as can be seen in Fig. 6.2 (??1 = 36 is similarly
chaotic for all ?).
0 2 4 6 8
W
5
10
15
<<
v
>>
n-1=36
n-1=18
n-1=6.3
Figure 6.3: The average velocity hhvii with varying ?. The error bars do not indicate
uncertainty in hhvii, they extend up and down by the amount hhvii .
85
0 2 4 6 8
W
50
100
150
200
250
<<
v
 
 
v
>>
n-1=36
n-1=18
n-1=6.3
Figure 6.4: The average inertial force hhv ?rvii with varying ?. The error bars
do not indicate uncertainty in hhv?rvii, they extend up and down by the amount
hhv?rvii .
0 2 4 6 8
W
100
200
300
400
<<
 
P
>>
n-1=36
n-1=18
n-1=6.3
Figure 6.5: The average pressure force hhrPii with varying ?. The error bars
do not indicate uncertainty in hhrPii, they extend up and down by the amount
hhrPii .
86
0 2 4 6 8
W
5
10
15
20
<<
v
 
F
>>
n-1=36
n-1=18
n-1=6.3
Figure 6.6: The average power hhv ? Fii with varying ?. The error bars do not
indicate uncertainty in hhv?Fii, they extend up and down by the amount hhv?Fii .
0 2 4 6 8
W
50
100
150
200
250
<<
H
>>
n-1=36
n-1=18
n-1=6.3
Figure 6.7: The average helicity hhHii = hhv?(r?v)ii with varying ?. The error
bars do not indicate uncertainty in hhHii, they extend up and down by the amount
ihHii .
87
0 2 4 6 8
W
2
4
6
8
<<
n
 
2
v
>>
n-1=36
n-1=18
n-1=6.3
Figure 6.8: The average viscous force hh?r2vii with varying ?. The error bars
do not indicate uncertainty in hh?r2vii, they extend up and down by the amount
hh?r2vii .
0 2 4 6 8
W
50
100
150
200
<<
2
W
?
v
>>
n-1=36n-1=18
n-1=6.3
Figure 6.9: The average Coriolis force hh2? ? vii with varying ?. The error bars
do not indicate uncertainty in hh2??vii, they extend up and down by the amount
hh2??vii .
88
If inertial waves are present in the most rapidly rotating  ows, we expect to
see frequencies near 2? in spectra of dynamical quantities. This signal is most easily
seen in the spectrum of the average power hv?Fi, shown in Fig. 6.10 for ? = 8. The
temporal Fourier transform of the average power is denoted Pow(!). The spectrum
consists of a peak near the highest allowed frequency (2? = 16) superimposed on a
background of turbulent noise. A Gaussian flt of this power spectrum to a function
of the form c + ae?b?(!?!?)2 to the range 11:6 < ! < 18:6 flnds the peak center
at !? = 14:3. Note that the k = 2^z + ^x and k = 2^z + ^y modes are expected to
have a frequency ! = 2p5? = 14:31. From equation (4.18), the expected quality
factor for k = 2^z + ^y and k = 2^z + ^x waves is Q = 25:8. However, a Gaussian
flt to the peak in Fig. 6.10 gives an inverse fractional full width at half height of
!=?! ? 6:5, approximately one fourth of the quality factor expected from linear
theory for perturbations about v = 0. The peak in Fig. 6.10 is a clear sign that
inertial waves are present.
89
5 10 15 20 25
w
10
20
30
40
50
60
Pow
H
w
L
Figure 6.10: The spectrum of the spatially averaged power hv ? Fi at the rotation
rate ? = 8 (and ??1 = 18). The peak is centered at ! = 14:3.
6.1.1 Flow Structure
The spatial structure of the  ow develops anisotropy as rotation rates increase
past ? = 2. Fig. 6.11 shows the structure of the pressure on the cube surface. The
trend toward two-dimensionalization of the  uid state is evident, as predicted by the
Taylor-Proudman Theorem. Figures 6.12 and 6.13 show the structure of the velocity
magnitude and Coriolis force magnitude on the cube surface. The  ow structure
undergoes a qualitative transition as the rotation rate increases past ? = 2.
90
Figure 6.11: The pressure P at the time t = T is shown on the cube surface for
various rotation rates and ??1 = 18. The scale has been normalized so that the color
spectrum approximately spans the range between the lowest and highest pressures
attained on the cube surface. For each state, the pressure increases from blue to
green to red, with pure green mapping to P = 0. The spatial average of the pressure
is deflned to be zero.
91
Figure 6.12: The velocity magnitude jvj at the time t = T is shown on the cube
surface for various rotation rates and ??1 = 18. The state is the same as that
shown in Fig. 6.11. The scale has been normalized so that the color spectrum
approximately spans the range between the lowest and highest velocity magnitudes
on the cube surface. The velocity increases from blue to green to red.
92
Figure 6.13: The magnitude of the Coriolis force j2??vj = 2?pv2x +v2y at the time
t = T is shown on the cube surface for various rotation rates and ??1 = 18. The
state is the same as that shown in Fig. 6.11 and 6.12. The scale has been normalized
so that the color spectrum spans the range between the lowest and highest Coriolis
force magnitude on the cube surface. The strength of the Coriolis force increases
from blue to green to red.
93
Figure 6.14: The velocity magnitude jvj at the time t = T is shown on the cube
surface for the rotation rates ? = 2 and ? = 4 for the high Reynolds number case
??1 = 36 at a high resolution (N=64). These similar-looking  ows are vastly difier-
ent in terms of their ability to generate a dynamo (as will be shown in Section 6.2).
The scale has been normalized so that the color spectrum approximately spans the
range between the lowest and highest velocity magnitudes on the cube surface. The
velocity increases from black to white.
6.1.2 The Energy Spectrum
The one dimensional time averaged kinetic energy spectrum is
Ev(k) = 1T ?t
?
Z T
t?
hE(k)v ijkj=kdt; (6.4)
where E(k)v = F
h
jvj2=2
i
, and h?ijkj=k indicates an average over all k vectors having
length k. As a test case we compute the energy spectrum for a high resolution test
case (N = 96, K = 48) with low viscosity ??1 = 50. Fig. 6.15 shows the results,
94
verifying that the code has produced a result roughly consistent with the expected
Kolmogorov result of Ev(k) ? k?5=3 in the inertial range. We regard the inertial
range as extending from k = 2 to k ? 10. Ev is largest for k = 1, which is the mode
that F injects energy into.
1 2 5 10 20 50
k
0.0001
0.001
0.01
0.1
1
E
v
H
k
L
 
E
v
-5 3
Figure 6.15: The normalized kinetic energy spectrum Ev(k)=Ev is shown for the test
case ??1 = 50, ? = 0, with resolution N = 96, ?t = 5 ? 10?4. The inertial range is
seen to follow the expected Kolmogorov scaling in which log
?
Ev(k)
?
?C??53log(k)
for some constant C?.
For the range of ? investigated, there is no substantial inertial range, the ar-
guments that lead to the k?2 spectrum in rotating turbulence do not hold here; a
much lower ? (requiring a much higher resolution) would be needed to see this be-
havior. Instead, the energy spectrum exhibits k?3 behavior as seen in [55], although
the conditions in this case are quite difierent (in [55] random forcing was applied at
small scales, and the k?3 behavior was seen for scales larger than the forcing scale).
95
Fig. 6.16 shows the efiect of rotation on the compensated kinetic energy spectrum
Ev(k)k3=Ev. These spectra are computed from K = 32 simulations.
2 5 10 20
k
0.15
0.2
0.3
0.5
0.7
1
1.5
E
v
H
k
L
k
3
 
E
v
W=0
W=1
W=2
W=4
W=6W=8
Figure 6.16: The compensated energy spectrum is shown for the rotation rates:
? = 0; 0:5; 1; 2; 4; 6; 8 with ??1 = 18.
One can identify a transition in the spectrum occurring between ? = 1 and
? = 6. Rotation rates below ? = 2 (? = 0; 0:5; 1) largely follow the non-rotating
expectation, with viscous dissipation at small scales causing the spectrum to curve
downward, deviating from a power law. Rotation rates above ? = 2 (? = 4; 6; 8)
approximately follow a k?3 law from k = 3 to k = 10. This transition in the
spectrum near ? = 2 coincides with the transition seen in the  ow structure images
of the previous section. High rotation rates have signiflcantly more energy in the
large scale k = 2 mode. We will see that the kinematic dynamics of the magnetic
fleld also change substantially as the rotation rate increases past ? = 2.
96
6.2 Kinematic Dynamo Characterization
The growth or decay of the magnetic energy is measured by the flnite time
exponential growth rate
 ? = 1?log
?E
B(t? +?)
EB(t?)
?
: (6.5)
The initial seed magnetic fleld B? at time t = t? is random, with a small energy
EB(t?) distributed evenly in the lowest 10 modes (k = 1 to k = 10). The initial
velocity fleld v? is taken from a long-time simulation and thus re ects a randomly
chosen state from the attractor of the rotating Navier-Stokes system (4.13). The
growth rate  ? appears to limit to a constant for large ?, we denote this constant
 1 as the average of  ? over random initial conditions for M long, (? = 300)
simulations. We have
 1 ? ? ? = 1M
MX
i=1
 i?; (6.6)
where  i? denotes the value of  ? for initial condition i (vi?, Bi?). The error in this
approximation due to flnite sampling (flnite M) is estimated as
? ? ? 1pM
vu
ut 1
M ?1
MX
i=1
( i? ? ^ ?)2; (6.7)
the sample standard deviation divided by pM. For ??1 = 18 we average over
M = 10 initial conditions, while for ??1 = 6:3 we typically average over more
(between M = 10 and M = 20). In the ??1 = 36 case, the resolution requirements
are such that we only average over M = 2 or 3 initial conditions. Figure 6.17 shows
a set of  i? time series; 10 traces for each of two difierent ? values. Dynamo action
is characterized by a positive  1 indicating average growth of the seed magnetic
97
fleld, while negative  1 indicates average decay of the fleld. Thus, in Figure 6.17
the black curves (??1 = 14) show dynamo action, while the red curves (??1 = 10)
show an absence of dynamo action.
50 100 150 200 250 300
t
-0.2
-0.1
0
0.1
0.2
s
t
Figure 6.17: An example of the magnetic energy growth rate  ? converging on  1
for several initial conditions. The red curves correspond to ??1 = 10 while the black
correspond to ??1 = 14. In all cases ? = 1, ??1 = 18. From these ensembles
we obtain the estimates  1 ? ? 300 ? ? 300 = ?0:0587 ? 0:0129 (??1 = 10) and
 1 ? ? 300 ?? 300 = 0:0967?0:0126 (??1 = 14).
For flxed ? and ?,  1 increases with ??1. Fig. 6.18 shows this typically linear
dependence. The critical value of ? for which  1 = 0 is denoted ?c. The  uid  ow
is independent of ?, and thus for flxed ? and ? the average velocity Urms = hhvii
is independent of ?. We may then deflne (again, for flxed ? and ?) the critical
magnetic Reynolds number as Rcm = 2?Urms=?c = 2?hhvii??1c .
98
11 11.5 12 12.5 13 13.5 14
h-1
-0.05
-0.025
0
0.025
0.05
0.075
0.1
s
300
Figure 6.18: The magnetic energy growth rate versus ??1 for ? = 18?1, ? = 0:5.
The red line is a linear flt to the data points, which intersects  300 = 0 at a critical
value of ?c ? 11:6?1.
To estimate ?c we flt the  1 values to a straight line (? 300 = a??1 +b). Giving
??1c = ?b=a. This flt takes into account the uncertainty ? 300 of each point by
flnding the line with maximum likelihood and assuming Gaussian statistics. This
procedure is performed for several values of ? and ?; the results are shown in
Fig. 6.19. Just as the energy spectrum and  ow structure change dramatically for
rotation rates above ? = 2, so does the propensity for dynamo action. The value of
??1c decreases dramatically as rotation increases, indicating that rotation is desirable
for the generation of a dynamo  ow. The observed transition in the kinetic energy
spectrum coincides with this drop in ??1c near ? = 2.
As rotation increases, so does hhvii, increasing the average input power hhv?Fii
99
proportionately. In an experiment the maximum input power is typically limited by
practical constraints. Thus, although the critical magnetic difiusivity has changed
favorably for higher rotation rates, the critical input power dependence on ? is cru-
cial. The critical magnetic Reynolds number Rcm is a more fair comparison because
it multiplies ??1c by the RMS velocity, which is roughly proportional to the power
input (see Fig. 6.6 and 6.3). Fig. 6.20 shows the critical magnetic Reynolds number
versus rotation rate. Although the drop in Rcm is not as large as for ??1c , there is
still a signiflcant drop as ? increases past ? = 2. Fig. 6.20 summarizes the main
result of this part of the thesis.
2 4 6 8
W
2
4
6
8
10
12
14
16
h
c-
1
n-1=36
n-1=18
n-1=6.3
Figure 6.19: The critical magnetic difiusivity ??1c as a function of rotation rate. The
error due to flnite sampling of initial conditions is on the order of the line thickness.
Faster rotation rates are seen to be conducive to dynamo action, with over a ten-fold
decrease in ??1c when ? > 4 compared to the no rotation case.
100
2 4 6 8
W
100
200
300
400
500
R
mc
n=36-1
n=18-1
n=6.3-1
Figure 6.20: The critical magnetic Reynolds number Rcm as a function of rotation
rate. Faster rotation rates are seen to be conducive to dynamo action. The flnite
sample size induced uncertainty is on the order of the line thickness.
The drop ofi of Rcm occurs for approximately the same values of ? for all
viscosities (??1 = 36, ??1 = 18, and ??1 = 6:3). Thus, the viscous force is unlikely
to be important in the observed behavior. The Rossby number, which compares the
strength of the Coriolis force to the strength of the inertial force (not the viscosity
as in the Ekman number), is the important dimensionless number. The favorable
conditions for dynamo action occur for ? > 4 corresponding approximately to Ro<
0:3 (for ??1 = 36 and ??1 = 18).
The spatial structure of the growing magnetic fleld is shown in the surface of
the cube for various rotation rates in Fig. 6.21. In these images, the values of ??1
are slightly above the value ??1c for the given rotation rate so that the magnetic
fleld experiences average growth. These images give a sense of the spatial scale
101
and distribution of growing of magnetic structures in each case. The qualitative
difierence for difierent ? largely re ects difierences due to the changing ?. Small
scale structure is heavily suppressed for smaller ??1.
Figure 6.21: The magnitude of the magnetic fleld jBj at the time t = T is shown
on the cube surface for various rotation rates. Each state occurs for a value of
? near (but below) ?c for that particular value of ?; from left to right they are
??1 = [12;14;14;9;2:5;1:2;0:8]. The scale has been normalized so that the color
spectrum spans the range between the lowest and highest value of jBj on the cube
surface, where the strength of the magnetic fleld increases from blue to green to red.
102
6.3 Connection to an Envisioned Experimental Situation
Although the studied system is highly idealized (periodic boundaries, external
ABC forcing, etc.), we are interested in how the reported results of section 6.2
can be used to inform and improve dynamo experiments. The important question
we wish to address is this: given an available average input power density ?P =
hh?F??vii, viscosity ??, and magnetic difiusivity ?? (and hence flxed Prm), is it generally
beneflcial to rotate the experimental apparatus (the goal being the generation of a
dynamo  ow)? In Fig. 6.20, the curves are for flxed ?, not flxed ?? (the relationship
between ? and ?? involves ?U and thus the forcing strength ?F), making Fig. 6.20
di?cult to use to answer the question posed above. In our envisioned experiment,
the power density ?P is controlled by the experimenter, who we suppose increases ?P
until a dynamo is achieved. Thus, we wish to locate the critical power density ?Pc,
above which a dynamo is obtained for various flxed rotation rates ??.
The average power density is
?P = hh?F? ?vii; (6.8)
which can also be expressed as
?P = ?F ?Up
3hhF?vii
= 2???
?U3
?L hhF?vii
= 2????L
 2???
?L?
?3
hhF?vii; (6.9)
where ?U depends on the magnetic difiusivity via ?U = (2???)=(?L?). Deflning the
103
dimensionless power to be P = ?P(?L=(2?))4(1=(????3)) we have
P = 1?3hhF?vii: (6.10)
This normalization is motivated by envisioning a situation in which a particular  uid
is used (hence ?? and ?? are flxed and known) in an experiment of flxed conflguration
(for us with our periodic boundary conditions this corresponds to having a flxed
known ?L). Thus, in such a case, the dimensional power density ?P is easily computed
from P, since ?L, ??, ??, and ?? are flxed and known.
From our normalized Navier-Stokes equation 4.12, we see that the average
value of hhF ? vii is in general a function of the dimensionless pair (?;?). The
normalized quantities ? and ? involve the dimensional variable ?U. We regard this as
undesirable since ?U depends on ?F which is not flxed in our envisioned experimental
situation. To remedy this, a new dimensionless variable is used,
C =
???L2
2??? =
2??
? ; (6.11)
which is the ratio of the magnetic difiusion time to the rotation period. Additionally,
Prm depends only on the  uid used; hence we express the power P as a function
of C, Prm, and ? to construct the space of dimensionless parameters (Prm;C;P)
which has three desirable properties: 1) this triple is uniquely determined from the
previously used dimensionless triple (?;?;?), 2) all three quantities are easily com-
puted from known dimensional quantities in an experiment, and 3) the relationship
between (??;??; ??) and (Prm;C;P) involves only flxed quantities.
Our procedure consists of flxing Prm and C, varying ? to numerically flnd
?c, and computing the critical power Pc(Prm;C) = P(Prm;C;?c). To accomplish
104
this using our results from the previous section and without performing additional
simulations for flxed values of C, we make the approximation that the surface sep-
arating the parameter space (?;?;?) into dynamo and non-dynamo regions can be
formed by linearly interpolating between the curves of Fig. 6.19. Fig. 6.22 shows
this surface with the addition of points at ??1 = 3. This approximation is justifled
by the observation that the curves in Fig. 6.19, being roughly similar, indicate a
smooth dependence on ?.
0
2
4
6
8
W
10
20
30
n-1
0
5
10
15
20
h-1
Figure 6.22: Surface seperating dynamo behavior from non-dynamo behavior for
the parameter region 3 ???1 ? 36, 0 ? ? ? 8. Points above this surface generate
a dynamo, while points below do not.
P is computed by flrst linearly interpolating between the curves of Fig. 6.6 to
obtain a value for hhF?vii for any given (?; ?) in the range 3 ???1 ? 36, 0 ? ? ? 8.
105
This value is then multiplied by ??3 to obtain P (of equation (6.10)). The surface in
Fig. 6.22 may then be mapped to the space (Prm;C;P). Fig. 6.23 shows an example
of how the critical surface maps to constant Prm planes, demonstrating that for a
flxed P, increasing the dimensionless rotation rate C is desirable for the generation
of a dynamo (in the region of Prm investigated here).
0 25 50 75 100 125
C
5
10
15
P
H
?
10
3
L
HC*,P*L
HC**,P**L
Dynamo
No dynamo
PcH--LPc
H-L
PcH+L
Figure 6.23: Curve separating dynamo behavior from non dynamo behavior in the
plane (P;C) for Prm = 0:5.
The top portion of the curve in Fig. 6.23 is denoted P(+)c , the critical power
above which a dynamo is always attained. The bottom surface is denoted P(?)c , it
represents the boundary of a small dynamo window which exists for low powers. In
addition, there is another critical power P(??)c <P(?)c such that no dynamo action
occurs for P <P(??)c . We deflne the point (C?;P?) where P(+)c and P(?)c meet; as
well as the point (C??;P??) where P(?)c and P(??)c meet (see Fig. 6.24). The curve
P(??)c is positive for all C >C?? (for zero power P we cannot have a dynamo). The
106
dynamo region is P >P(+)c when C <C??, P(??)c >P >P(?)c and P >P(+)c when
C?? <C <C?, and P >P(??)c when C >C?.
The dynamo window that exists in the range P(??)c > P > P(?)c is shown in
Fig. 6.24 along with the curve PTc for which the  ow has complex time dependence
for P >PTc and is time-independent or periodic for P <PTc . The curve P(??)c can
be seen to exist in a region in which the power is so low that the resulting  ow is not
turbulent. The non-turbulent dynamo region evidently yields time independent or
periodic  ows having a spatial structure particularly e?cient for dynamo generation.
For increasing powers, turbulence develops above PTc in the form of  uctuations on
top of a mean  ow. This mean  ow is similar in spatial structure to the time-
independent  ows just below PTc . Our interpretation of Fig. 6.24 is that when
the additional turbulence is weak, the mean  ow is able to preserve the e?cient
dynamo; however, as P increases past P(?)C , turbulence causes the  ow to become
disordered (perhaps introducing a greatly enhanced efiective turbulence induced
magnetic difiusivity) and the dynamo is lost. Above P(+)c the power is su?ciently
high that dynamo action occurs for highly turbulent  ow.
107
44 48 52 56 60C
25
50
75
100
125
P
Dynamo
No dynamo
No dynamo
PcT
PcH-L
PcH--L
HC**,P**L
Figure 6.24: The critical power curves P(??)c , P(?)c , and PTc for Prm = 0:2.
Experiments operate in the regime Prm ? 1. In principle, we wish to locate
the critical power curves P(+)c , P(?)c , and P(??)c in the limit as Prm goes to zero.
The increase in ??1c with ??1 seen in Fig. 6.19 corresponds to a rising P(+)c as Prm is
decreased. As one approaches the low Prm limit, we expect that P(+)c will limit to a
curve denoted P(+)c1 corresponding to a ?-independent ?c. This expectation is based
on results obtained in the non-rotating case, demonstrating that Rcm increases for
decreasing Prm, reaching a plateau for small Prm (Prm ? 0:3) [40]. We expect that
the shape of P(+)c will not be altered drastically in the limit Prm ! 0. However,
investigation of very small magnetic Prandtl numbers proved to be beyond the
capabilities of our available computer resources. The lowest Prm for which we have
obtained the entire upper curve P(+)c is Prm ? 0:45. (data for which ??1 > 36 would
be required to obtain P(+)c for lower Prm).
108
Although we only obtain the upper curve for a relatively small region of Prm,
we are able to obtain large portions of the P(?)c curve as low as Prm ? 0:1. We are
uncertain as to the behavior of P(?)c in the limit Prm ! 0; its behavior as Prm is
decreased is shown in Fig. 6.25 (the dynamo window bounded by P(?)c and P(??)C
exists for Prm as low as Prm = 0:1 and thus we cannot rule out its existence for
Prm ! 0). The curve P(??)c versus C is relatively unchanged as Prm decreases.
0 25 50 75 100 125
C
5
10
15
P
H
?
10
3
L
Prm= 0.8
Prm= 0.6
Prm= 0.5
Prm= 0.4 Prm= 0.3
Prm= 0.2
Figure 6.25: Critical power curves P(+)c and P(?)c versus C for several values of Prm;
they are Prm = 0:8; 0:6; 0:5; 0:4; 0:3; 0:2. Curves P(+)c which do not extend back
to C = 0, are incomplete because it would require data from outside of the region
3 ???1 ? 36.
109
6.4 Conclusions
As early as 1970 it was hypothesized that rotation increases the propensity
for dynamo action in a magnetohydrodynamic  uid [53, 54]. Here we have shown,
through the use of numerical simulations, that rotation is desirable for dynamo
action. Moreover, the behavior change is marked by a sharp transition in the critical
magnetic Reynolds number. This transition in magnetic behavior was found to
coincide with changes in the  ow structure and kinetic energy spectrum. In addition,
the  ow helicity grows with rotation rate, providing a possible explanation for the
observed magnetic behavior.
6.4.1 Future Prospects
Our study has focused on external forcing on the scale of the domain (the
largest scale k = 1 mode). However, quasi two-dimensional  ows, which occur for
the highest rotation rates investigated, are known to have an inverse energy cascade
to larger scales. Studying rotation with a forcing scale small compared to the domain
(perhaps with the use of a Taylor-Green vortex with k? > 1) is thus of great interest.
It is also of interest how boundaries efiect the  ow. These topics are appropriate
for future study, when more computational recourses become available.
110
Appendix A
Geometric Optics Derivation of the Shadowgraph Light Intensity
The path of a light ray parameterized by the path length s is denoted r(s).
Using this parameterization, v ? dr=ds has magnitude 1. A light ray beginning
at r? = x?^x + y?^y in the xy plane with initial v vector pointing down along the
z-axis (v? = ?^z), travels into a  uid contained between the planes z = d=2 and
z = ?d=2. Upon hitting the z = ?d=2 plane, the light ray re ects and exits the
 uid for subsequent imaging (through a series of lenses followed by a CCD).
In a medium with general index of refraction n(x;y;z) the path satisfles
d
ds
?
n(r(s))dr(s)ds
?
= rn(r(s)): (A.1)
Applying the chain rule we have
 
rn(r)? drds
? dr
ds +n(r)
d2r
d2s = rn(r): (A.2)
Substituting v = dr=ds leads to
dv
ds = R(r)?(R(r)?v)v; (A.3)
dr
ds = v; (A.4)
where
R(r) = rn(r)n(r) : (A.5)
The right hand side of equation (A.3) is simply the projection of R(r) onto the
plane perpendicular to v. The index of refraction of a  uid can be expressed as
111
n(x;y;z) = 1 +?n(x;y;z) where ?n is treated as a small perturbation, ?n?n. We
wish to treat the problem to flrst order in ?n. Thus R(r) becomes
R(r) ?r?n(r): (A.6)
To lowest order we ignore ray deviations within the  uid so that the ray exits
the  uid at the same point it entered (r?). Additionally, the vector v is only slightly
modifled during transit through the  uid, thus we substitute v = ?^z in the right
hand side of (A.3). This gives the simplifled system
dv
ds = R?(r): (A.7)
Here, the subscript ? indicates the projection of the R vector onto the xy plane.
The exiting v vector ve = ?vx^x + ?vy ^y + ^z picks up a small component in the
xy plane. The ray then travels in a straight line along the direction ve to a flnal
location r = (x? +?x)^x + (y? +?y)^y + (z1 +d=2)^z in the image plane, a distance
z1  d above the  uid. The ofisets ?x = z1?vx, and ?y = z1?vy are treated as
perturbations. Integrating (A.7) gives
?v? ? ?vx^x +?vy ^y
=
Z ?d=2
d=2
R?(r?)ds+
Z d=2
?d=2
R?(r?)ds
= 2r?
Z d=2
?d=2
?n(r? +z^z)dz (A.8)
?r? ? ?x^x +?y^y = z1?v?
= 2z1r?
Z d=2
?d=2
?n(r? +z^z)dz
= 2z1dr??n(r?); (A.9)
112
where the overbar indicates vertically averaging, ?n = (1=d)Rd=2?d=2?ndz. The map
from the ray entrance point r? to the flnal horizontal location rf where it hits the
image plane is
rf = P(r?) ? r? + 2z1dr??n(r?): (A.10)
The map P is one-to-one since the deviation ?r? is small. The determinant of the
Jacobian is the factor by which areas grow under the action of the map P:
jDP(r?)j = drfxdr
?x
drfy
dr?y ?
drfx
dr?y
drfy
dr?x
=
 
1 + 2z1dd
2?n
d2x (r?)
? 
1 + 2z1dd
2?n
d2y (r?)
?
?
4z21d2
 d2?n
dxdy(r?)
?2
? 1 + 2z1dr2??n(r?): (A.11)
Equation (A.11) is accurate to flrst order in ?n. For an incident light intensity I?
at the point r? the point rf will receive intensity I = jDP(r?)j?1I?. Thus we have
the relation
I(x;y) = I?(x;y)1 + 2z
1dr2??n(x;y)
: (A.12)
The  uid?s index of refraction varies due to spatial temperature variations.
For slight dependence of n on the temperature deviation  (ignoring the higher
derivatives such as d2n=dT2) we have: ?n = ?n? +(dn=dT) . Typically dn=dT < 0;
a convenient form for the flnal shadowgraph formula is then
I(x;y) = I?(x;y)1?2z
1djdn=dTjr2?? (x;y)
: (A.13)
113
Appendix B
Short-Time Nonlinear Evolution of an Initially Gaussian PDF
The goal of this appendix is to compute the initial time derivative of the mean,
covariance, skewness, and kurtosis of a multivariate normal distribution (PDF, or
equivalently, density of an ensemble of states) under nonlinear evolution. The last
two quantities (skewness and kurtosis) are measures of deviations from normality.
We consider a continuous time dynamical system
@x
@t = F(x); (B.1)
in which the N components of F can be written
Fi(x) = Ci +Aijxj +Bijkxjxk; (B.2)
where repeated indices are summed. Without loss of generality, the B tensor is sym-
metric in its last two indices; Bijk = Bikj. This form of F(x), in which second order
terms are the highest present, is very common. The Boussinesq equations (1.1), the
standard Navier-Stokes equations (4.12), and the MHD equations (4.1) and (4.2)
take this form when discretized in space, with all nonlinearities arising from terms
which are second order in their state variables.
Denote the PDF as ?(x). Conservation of probability leads to the continuity
equation in state space
@?
@t +r?(?F) =
@?
@t +?r?F +r??F = 0: (B.3)
114
The initial (t = t?) density ??(x) is Gaussian with covariance ?? and mean ??,
??(x) = 1(2?)N=2j?
?j1=2
exp
?
?12(x???)T??1? (x???)
?
: (B.4)
For this distribution we have r?? = ???1? (x ? ??)??. It will be convenient to
transform into a coordinate system in which ?? is centered on the origin and nor-
mally distributed with covariance matrix equal to the N?N identity matrix. This
coordinate transform is given by
? = P?1? (x???); (B.5)
where P is real, symmetric, and satisfles
P2? = ??: (B.6)
Using these new deflnitions, r?? = ?P?1? ???. The initial time evolution of the
density is then
_?? ? @?@t
flfl
flfl
t=t?
= P?1? ??? ?F???r?F: (B.7)
Using (B.4) and (B.2), equation (B.7) becomes
_?? = ???(Ci +Aijxj +Bijkxjxk)P?1?il ?l ?Aii ?2Biijxj?: (B.8)
Eliminating x using the substitution x = ?? + P?? leads to the form
_?? = ??
h
CiP?1?il ?l +Aij??jP?1?il ?l +AijP?jpP?1?il ?p?l +
Bijk??j??kP?1?il ?l + 2Bijk??kP?jpP?1?il ?p?l +BijkP?knP?jpP?1?il ?n?p?l ?
Aii ?2Biij??j ?2BiijP?jl?l
i
: (B.9)
115
Although long and seemingly cumbersome, this form is convenient because we can
easily make use of the relations
Z
??dx = 1 ; (B.10)
Z
?i??dx = 0 ; (B.11)
Z
?i?j??dx = ?ij ; (B.12)
Z
?i?j?k??dx = 0 ; (B.13)
Z
?i?j?k?l??dx = ?ij?kl +?ik?jl +?il?jk ;
? Wijkl f??g ; (B.14)
Z
?i?j?k?l?m??dx = 0 ; (B.15)
Z
?i?j?k?l?m?n??dx = ?ijWklmn f??g+?ikWjlmn f??g+
?ilWjkmn f??g+?imWjkln f??g+?inWjklm f??g ;
? Wijklmn f???g ; (B.16)
Z
?i?j?k?l?m?n?p??dx = 0 ; (B.17)
...
These relations follow from noting that j??j1=2 = jP?j and under change of variables
dx !jP?jd?.
At any time, the distribution mean is computed as
? =
Z
x?(x)dx; (B.18)
and its initial time derivative is
_?? ? @?@t
flfl
flfl
t=t?
=
Z
x _??(x)dx =
Z
?? _??(x)dx +
Z
P?? _??(x)dx: (B.19)
116
Noting that R _??(x)dx = 0 we have
_??i = P?ij
Z
?j _??dx: (B.20)
Only terms with odd powers of ? in (B.9) will give non-zero contributions to this
integral. After plugging (B.9) into (B.20) and using the relations P?ij = P?ji, ??ij =
P?ikP?jk, and P?ijP?1?jk = ?ik, we obtain
_?i = Fi(?) + ??jkBijk: (B.21)
Evidently the mean evolves initially just as if it were a state evolving under the action
of F, except for a nonlinear correction due to the flnite extent of the distribution.
For very tight distributions (small ??), the distribution mean follows _? = F(?).
The distribution covariance is given by
?ij =
Z
(xi ??i)(xj ??j)?(x)dx; (B.22)
and its initial time derivative is
_??ij ? @?ij
@t
flfl
flfl
t=t?
=
Z
(xi ???i)(xj ???j) _??dx?
_??i
Z
(xj ???j)??dx? _??j
Z
(xi ???i)??dx
= P?inP?jm
Z
?n?m _??dx?( _??iP?jn + _??jP?in)
Z
?n??dx
= P?inP?jm
Z
?n?m _??dx: (B.23)
Only terms with even powers of ? in (B.9) will give non-zero contributions to this
integral. After plugging (B.9) into (B.23) we obtain
_?? = DF??? +?DF????T; (B.24)
117
where DF?ij ?Aij + 2Bijk??k is deflned as the Jacobian of F at the point ??. The
same formula (B.24) is obtained when the dynamics are linear.
The skewness and kurtosis measures are given by [60] as
s =
Z ?
(x??)T??1(y ??)?3?(x)?(y)dxdy; (B.25)
? =
Z ?
(x??)T??1(x??)?2?(x)dx; (B.26)
Note that the skewness of ?? is zero and the kurtosis of ?? is N2 + 2N (recall that
N is the dimension of the space).
The initial time derivative of the skewness is
_s? ? @s@t
flfl
flfl
t=t?
=
Z ?
(x???)T??1? (y ???)?3 ( _??(x)??(y) +??(x) _??(y))dxdy ?
3
Z ?
(x???)T??1? (y ???)?2
?
(x???)T??1? _????1? (y ???) +
_?T? ??1? (y ???) + _?T? ??1? (x???)
?
??(x)??(y)dxdy: (B.27)
Where we have used the relation
@
@t?
?1
flfl
flfl
t=t?
= ???1? _????1? : (B.28)
The skewness is symmetric with respect to the interchange of x and y. Using this
fact, along with the notation ?(x) = P?1? (x???), ?(y) = P?1? (y???), and deflning
L ? P?1? DF?P? + P?DF T? P?1? ; (B.29)
we get
_s? = 2
Z h
(?(x))T?(y)
i3
_??(x)??(y)dxdy ?
3
Z h
(?(x))T?(y)
i2
(?(x))T L?(y)??(x)??(y)dxdy: (B.30)
118
Both integrals have odd powers of ?y, thus evaluating the dy portion of both
integrals flrst gives _s? = 0.
The initial time derivative of the kurtosis is
_?? ? @?@t
flfl
flfl
t=t?
=
Z ?
(x???)T??1? (x???)?2 _??dx?
2
Z ?
(x???)T??1? (x???)?
?
(x???)T??1? _????1? (x???) +
2 _?T? ??1? (x???)
?
=
Z ?
?T??2 _??dx?2
Z ?
?T???T L???dx
=
Z
?i?i?j?j _??dx?2
Z
Ljl?i?i?j?l??dx: (B.31)
Using (B.9) and the fact that ?ii = N, we arrive at
_?? = 4(N + 2)DF?ii ?2(N + 2)Lii
= 0;
since Lii = 2DF?ii (which can be shown using (B.29)).
In summary, we computed the short-time evolution of the mean, covariance,
skewness, and kurtosis for a multivariate normal distribution. The results are sum-
marized as
_? = F(?) + B??; (B.32)
_?? = DF??? +?DF????T; (B.33)
_s? = 0; (B.34)
_?? = 0: (B.35)
The initial growth of these measures of skewness and kurtosis are zero; thus the
119
growth of skewness and kurtosis must be (at least) second order in time for short
times.
120
Bibliography
[1] F. T. Arecchi, S. Boccaletti, and P. Ramazza. Pattern formation and competi-
tion in nonlinear optics. Physics Reports, 318:1{2, September 1999.
[2] H. L. Swinney and V. I. Krinsky. Waves and Patterns in Chemical and Biolog-
ical Media. The MIT Press, December 1991.
[3] Stephen W. Morris, Eberhard Bodenschatz, David S. Cannell, and Guenter
Ahlers. Spiral defect chaos in large aspect ratio rayleigh-b?enard convection.
Phys. Rev. Lett., 71(13):2026{2029, Sep 1993.
[4] M. C. Cross and P. C. Hohenberg. Pattern formation outside of equilibrium.
Rev. Mod. Phys., 65(3):851{1086, 1993.
[5] P. C. Hohenberg and B. I. Shraiman. Chaotic behavior of an extended system.
Physica D Nonlinear Phenomena, 37:109{115, July 1989.
[6] Hao-wen Xi, J. D. Gunton, and Jorge Vi~nals. Spiral defect chaos in a model of
rayleigh-b?enard convection. Phys. Rev. Lett., 71(13):2030{2033, Sep 1993.
[7] R.E. Kalman. A new approach to linear flltering and prediction problems. J.
Basic Eng., 82:35{45, 1960.
[8] Dan Simon. Optimal State Estimation: Kalman, H Inflnity, and Nonlinear
Approaches. Wiley-Interscience, 2006.
[9] Geir Evensen. Data Assimilation: The Ensemble Kalman Filter. Springer,
2006.
[10] E. Ott, B. R. Hunt, I. Szunyogh, A. V. Zimin, E. J. Kostelich, M. Corazza,
E. Kalnay, D. J. Patil, and J. A. Yorke. Estimating the state of large spatio-
temporally chaotic systems. Physics Letters A, 330:365{370, September 2004.
[11] E. Ott, B. R. Hunt, I. Szunyogh, A. V. Zimin, E. J. Kostelich, M. Corazza,
E. Kalnay, D. J. Patil, and J. A. Yorke. A local ensemble Kalman fllter for
atmospheric data assimilation. Tellus Series A, 56:415{+, October 2004.
[12] B.R. Hunt, E.J. Kostelich, and I. Szunyogh. E?cient data assimilation for
spatiotemporal chaos: A local ensemble transform kalman fllter. Physica D,
230:112{126, 2007.
[13] Jefirey S. Whitaker and Thomas M. Hamill. Ensemble data assimilation with-
out perturbed observations. Mon. Wea. Rev., 130:1913, 2002.
[14] M. K. Tippett et al., Mon. Wea. Rev. 131? , 1485 (2003).
121
[15] Craig H. Bishop, Brian J. Etherton, and Sharanya J. Majumdar. Adaptive
sampling with the ensemble transform kalman fllter. part i: Theoretical aspects.
Mon. Wea. Rev., 129:420, 2001.
[16] J.L. Anderson. An ensemble adjustment kalman fllter for data assimilation.
Mon. Wea. Rev., 129:2884, 2001.
[17] J.S. Whitaker T.M. Hamill and C. Snyder. Distance-dependent flltering of
background error covariance estimates in an ensemble kalman fllter. Mon. Wea.
Rev., 129:2776{2790, 2001.
[18] P. L. Houtekamer and Herschel L. Mitchell. A sequential ensemble kalman fllter
for atmospheric data assimilation. Mon. Wea. Rev., 129:123, 2001.
[19] E. Bodenschatz, W. Pesch, and G. Ahlers. Recent Developments in Rayleigh-
B?enard Convection. Annual Review of Fluid Mechanics, 32:709{778, 2000.
[20] D. A. Egolf, I. V. Melnikov, W. Pesch, and R. E. Ecke. Mechanisms of exten-
sive spatiotemporal chaos in Rayleigh-B?enard convection. Nature, 404:733{736,
April 2000.
[21] K.-H. Chiam, M. C. Cross, H. S. Greenside, and P. F. Fischer. Enhanced tracer
transport by the spiral defect chaos state of a convecting  uid. Phys. Rev. E,
71(3):036205{+, March 2005.
[22] R. Schmitz, W. Pesch, and W. Zimmermann. Spiral-defect chaos: Swift-
Hohenberg model versus Boussinesq equations. Phys. Rev. E, 65(3):037302{+,
March 2002.
[23] M. C. Cross and Y. Tu. Defect Dynamics for Spiral Chaos in Rayleigh-B?enard
Convection. Physical Review Letters, 75:834{837, July 1995.
[24] K.-H. Chiam, M. R. Paul, M. C. Cross, and H. S. Greenside. Mean  ow and
spiral defect chaos in rayleigh-b?enard convection. Phys. Rev. E, 67(5):056206,
May 2003.
[25] Y. Hu, R. Ecke, and G. Ahlers. Convection for Prandtl numbers near 1: Dy-
namics of textured patterns. Phys. Rev. E, 51:3263{3279, April 1995.
[26] M. Paul and M. Einarsson. Extensive Chaos in Rayleigh-Benard Convection.
APS Meeting Abstracts, pages A10+, November 2006.
[27] F. H. Busse. Non-linear properties of thermal convection. Reports of Progress
in Physics, 41:1929{1967, December 1978.
[28] L. S. Tuckerman. Divergence-free velocity flelds in nonperiodic geometries. J.
Comput. Phys., 80(2):403{441, 1989.
[29] J. Swift and P. C. Hohenberg. Hydrodynamic  uctuations at the convective
instability. Phys. Rev. A, 15(1):319{328, Jan 1977.
122
[30] M. C. Cross and P. C. Hohenberg. Pattern formation outside of equilibrium.
Rev. Mod. Phys., 65(3):851, Jul 1993.
[31] J. Lega, J. V. Moloney, and A. C. Newell. Swift-hohenberg equation for lasers.
Phys. Rev. Lett., 73(22):2978{2981, Nov 1994.
[32] W. Merzkirch. Flow visualization. New York, Academic Press, Inc., 1974. 258
p., 1974.
[33] S. Rasenat, G. Hartung, B. L. Winkler, and I. Rehberg. The shadowgraph
method in convection experiments. Exp. Fluids, 7:412{420, 1989.
[34] J.L. Anderson and S.L. Anderson. A monte carlo implementation of the non-
linear flltering problem to produce ensemble assimilations and forecasts. Mon.
Wea. Rev., 127:2741, 1999.
[35] Louis M. Pecora and Thomas L. Carroll. Synchronization in chaotic systems.
Phys. Rev. Lett., 64(8):821{824, Feb 1990.
[36] M. R. Paul, M. I. Einarsson, P. F. Fischer, and M. C. Cross. Extensive chaos in
rayleigh-b[e-acute]nard convection. Physical Review E (Statistical, Nonlinear,
and Soft Matter Physics), 75(4):045203, 2007.
[37] S. Morris, E. Bodenschatz, D. Cannell, and G. Ahlers. The spatiotemporal
structure of spiral-defect chaos. Physica D, 97:164, 1996.
[38] M. Oczkowski, I. Szunyogh, and D. J. Patil. Mechanisms for the Development
of Locally Low-Dimensional Atmospheric Dynamics. Journal of Atmospheric
Sciences, 62:1135{1156, April 2005.
[39] Y. Ponty, H. Politano, and J.-F. Pinton. Simulation of Induction at Low Mag-
netic Prandtl Number. Physical Review Letters, 92(14):144503{+, April 2004.
[40] Y. Ponty, P. D. Mininni, D. C. Montgomery, J.-F. Pinton, H. Politano, and
A. Pouquet. Numerical Study of Dynamo Action at Low Magnetic Prandtl
Numbers. Physical Review Letters, 94(16):164502{+, April 2005.
[41] Y. Ponty, P. D. Mininni, J.-F. Pinton, H. Politano, and A. Pouquet. Dynamo
action at low magnetic Prandtl numbers: mean  ow versus fully turbulent
motions. New Journal of Physics, 9:296{+, August 2007.
[42] D. Sweet, E. Ott, J. M. Finn, T. M. Antonsen, and D. P. Lathrop. Blowout
bifurcations and the onset of magnetic activity in turbulent dynamos. Physical
Review E, 63(6):066211{+, June 2001.
[43] N. Seehafer, F. Feudel, and O. Schmidtmann. Nonlinear dynamo with ABC
forcing. Astronomy and Astrophysics, 314:693{699, October 1996.
[44] Paul H. Roberts and Gary A. Glatzmaier. Geodynamo theory and simulations.
Rev. Mod. Phys., 72(4):1081{1123, Oct 2000.
123
[45] C. A. Jones and P. H. Roberts. Convection-driven dynamos in a rotating plane
layer. Journal of Fluid Mechanics, 404:311{343, February 2000.
[46] M. Meneguzzi and A. Pouquet. Turbulent dynamos driven by convection. Jour-
nal of Fluid Mechanics, 205:297{318, 1989.
[47] A. Brandenburg, R.L. Jennings, ?A. Nordlund, M. Rieutord, R.F. Stein, and
I. Tuominen. Magnetic structures in a dynamo simulation. J. Fluid Mech.,
306:325{352, 1996.
[48] K. R. Lang. Astrophysical Data I. Planets and Stars. Astrophysical Data
I. Planets and Stars, X, 937 pp. 33 flgs.. Springer-Verlag Berlin Heidelberg
New York, 1992.
[49] P. A. Davidson. An Introduction to Magnetohydrodynamics. An Introduction to
Magnetohydrodynamics, by P. A. Davidson, pp. 452. ISBN 0521791499. Cam-
bridge, UK: Cambridge University Press, March 2001., March 2001.
[50] F. Krause and K.-H. Raedler. Mean-fleld magnetohydrodynamics and dynamo
theory. Oxford, Pergamon Press, Ltd., 1980. 271 p., 1980.
[51] D. Galloway and U. Frisch. A note on the stability of a family of space-periodic
Beltrami  ows. Journal of Fluid Mechanics, 180:557{564, 1987.
[52] H. K. Mofiatt. Magnetic fleld generation in electrically conducting  uids. Cam-
bridge, England, Cambridge University Press, 1978. 353 p., 1978.
[53] H. K. Mofiatt. An approach to a dynamic theory of dynamo action in a rotating
conducting  uid. Journal of Fluid Mechanics, 53:385{399, 1972.
[54] H. K. Mofiatt. Dynamo action associated with random inertial waves in a
rotating conducting  uid. Journal of Fluid Mechanics, 44:705{719, 1970.
[55] L. M. Smith and F. Walefie. Transfer of energy to two-dimensional large scales
in forced, rotating three-dimensional turbulence. Physics of Fluids, 11:1608{
1622, June 1999.
[56] P. K. Yeung and Y. Zhou. Numerical study of rotating turbulence with external
forcing. Physics of Fluids, 10:2895{2909, November 1998.
[57] J?er^ome Paret and Patrick Tabeling. Experimental observation of the two-
dimensional inverse energy cascade. Phys. Rev. Lett., 79(21):4162{4165, Nov
1997.
[58] M. Rieutord and L. Valdettaro. Inertial waves in a rotating spherical shell.
Journal of Fluid Mechanics, 341:77{99, June 1997.
[59] The library used for fast Fourier transforms is the fitw library
(http://www.fitw.org/).
124
[60] K. V. MARDIA. Measures of multivariate skewness and kurtosis with applica-
tions. Biometrika, 57(3):519{530, 1970.
125