# Shadowing the rotating annulus. Part I:
Measuring candidate trajectory shadowing times

Roland M. B. Young Roman Binter Falk Niehörster Centre for the Analysis of
Time Series, London School of Economics, London, UK  Atmospheric, Oceanic and
Planetary Physics, Department of Physics, University of Oxford, Oxford, UK
Bermuda Institute of Ocean Sciences, St. George’s, Bermuda

###### Abstract

An intuitively necessary requirement of models used to provide forecasts of a
system’s future is the existence of shadowing trajectories that are consistent
with past observations of the system: given a system-model pair, do model
trajectories exist that stay reasonably close to a sequence of observations of
the system? Techniques for finding such trajectories are well-understood in
low-dimensional systems, but there is significant interest in their
application to high-dimensional weather and climate models. We build on work
by Smith et al. [2010, Phys. Lett. A, 374, 2618–2623] and develop a method for
measuring the time that individual “candidate” trajectories of high-
dimensional models shadow observations, using a model of the thermally-driven
rotating annulus in the perfect model scenario. Models of the annulus are
intermediate in complexity between low-dimensional systems and global
atmospheric models. We demonstrate our method by measuring shadowing times
against artificially-generated observations for candidate trajectories
beginning a fixed distance from truth in one of the annulus’ chaotic flow
regimes. The distribution of candidate shadowing times we calculated using our
method corresponds closely to (1) the range of times over which the
trajectories visually diverge from the observations and (2) the divergence
time using a simple metric based on the distance between model trajectory and
observations. An empirical relationship between the expected candidate
shadowing times and the initial distance from truth confirms that the method
behaves reasonably as parameters are varied.

This paper was submitted to Physica D in 2010, but, after review, was not
accepted. We no longer have the time or resources to work on this topic, but
would like this record of our work to be available for others to read, cite,
and follow up.

###### keywords:

Shadowing trajectories; Rotating annulus; Dimensional scaling; Weather and
climate forecasting; Chaotic dynamics; Perfect Model Scenario.

††journal: arXiv

##  1 Introduction

In forecasting one is often forced to use complex numerical models to predict
the future behaviour of physical or mathematical systems. A good example are
General Circulation Models (GCMs), which are used to predict the future
behaviour of the atmosphere. GCMs have state spaces of over ten million
dimensions and integrate differential equations describing complex dynamical
and chaotic behaviour. Decisions on how to act over future timescales of days
(in the case of weather) up to decades (in the case of climate) are made based
on the output from these models, and so decision makers desire these forecasts
to be informative.

An intuitive requirement of models used to provide forecasts of a system’s
future is the existence of model trajectories that are consistent with past
observations of the system. Given a system-model pair, do model trajectories
exist that stay reasonably close to a sequence of observations of the system?
Such trajectories are said to shadow the observations. While the existence of
a shadowing trajectory is not a sufficient condition for successful forecasts,
it is intuitively a necessary one.

The concept of shadowing appears in several slightly different contexts, and
it is important to differentiate between these. In its original context, so-
called $\epsilon$-shadowing, the shadowing problem is concerned with whether
there exists a true solution to a differential equation that remains
sufficiently close to a given numerical solution. If no true solution exists
that shadows the numerical solution, can we learn anything at all about the
true solution from the numerical solution? The $\epsilon$-shadowing time is
formalised in the Anosov-Bowen shadowing lemma (Bowen, 1975). Second, a
superficially similar but fundamentally different context is $\iota$-shadowing
(Gilmour, 1998). This is concerned with whether a trajectory of a numerical
model exists that shadows a given sequence of observations of a system. The
$\iota$-shadowing time is the maximum model trajectory length that, among all
possible trajectories, remains consistent with the observational noise, and is
defined with respect to a particular noise model (Gilmour, 1998, p. 47). In
this context, if the model does not admit such a trajectory then there is no
such initial condition that remains close to the observations, and the problem
is then one of model error, not chaos. Third and finally, the $\phi$-shadowing
problem (Smith, 2000, p. 52) is concerned with whether there exists a model
trajectory that remains close enough to a given set of observations for
“useful” forecasts to be possible. As what is “useful” depends on context,
this use of the term is necessarily subjective. For example, one could argue
that a GCM solution that predicts the wrong amount of rain in a given place
but does predict that rain will fall is “useful” enough to constitute a
$\phi$-shadowing trajectory.

In the atmospheric context we are interested in how well particular models
shadow observations. Several methods exist that attempt to find shadowing
trajectories directly from a sequence of observations (e.g. Gilmour, 1998;
Stemler and Judd, 2009), and these methods are well-understood using low-
dimensional (low-$N$) systems such as the Lorenz equations (Lorenz, 1963)
($N=3$) and the Ikeda map (Ikeda, 1979) ($N=2$). Whether these methods are
feasible when applied to high-dimensional models of the atmosphere remains
poorly understood, however, as is the time over which atmospheric models can
shadow past observations.

Figure 1: Schematic of the thermally-driven rotating annulus experiment. The
inner cylinder is at radius $a=2.5$ cm, the outer cylinder is at radius
$b=8.0$ cm, and the cylinders have a depth of $d=14.0$ cm. The fluid between
the cylinders is a mixture of glycerol and water, 17%/83% by volume. The
cylinders are maintained at constant temperatures $T_{a}=18$ ∘C and $T_{b}=22$
∘C respectively, giving a constant temperature difference across the fluid of
$\Delta T=4$ degC. The whole system is rotated about a vertical axis with
angular velocity $\Omega=1.00$ rad s-1.

Intermediate in complexity between low-dimensional systems and the atmosphere
are laboratory experiments, such as the rotating annulus (Hide, 1953; Hide and
Mason, 1975; Hignett et al., 1985; Read et al., 1992). This system consists of
a fluid contained between two concentric cylinders, rotating about a vertical
axis with a horizontal temperature gradient between the cylinders (Fig. 1).
This setup mimics the major influences acting on a planet’s atmosphere: the
effects of rotation, gravity, and the temperature difference between low and
high latitudes. The rotating annulus has been well-established over 50 years
as a useful analogue for certain types of atmospheric dynamical behaviour,
especially the general circulation of the atmosphere at mid-latitudes (Winston
and Krueger, 1961; Lorenz, 1967). Models used to simulate the behaviour of
fluid in the annulus are intermediate in complexity between low-dimensional
systems and global atmospheric models; the dimension of these models are
$O(10^{4})$, compared with $O(10)$ for low-dimensional systems and $O(10^{7})$
for GCMs. Experimental systems are also intermediate in their physical
idealization of the system being studied, between GCMs with a low degree of
idealization, and systems such as the Lorenz equations, which describe a
highly-idealized representation of convective rolls at the onset of
convection. This paper is the first in a larger body of work exploring
shadowing techniques using the rotating annulus. We aim to understand
shadowing methods better by examining them in a real but idealised system, and
also improve our understanding of the annulus and the models used to study it
through their shadowing properties.

Later work will address the problem of finding the longest model trajectories
that shadow a given sequence of observations. Here we focus on how to measure
the shadowing time itself. Most methods that look for shadowing
trajectories111Hammel et al. (1988) is an exception, but would be very
difficult in practice to extend to a high-dimensional system. generate a set
of candidate trajectories, each of which may or may not shadow the
observations. One then has to measure how long each candidate shadows, and the
longest of these candidate shadowing times gives a lower bound on the true
shadowing time for that model — observation pair. The main aim of this paper
is to present a method for identifying shadowing times that can be used with
high-dimensional systems, like the rotating annulus, when the observational
error distribution is unbounded. We begin from a definition outlined in Smith
et al. (2010) which has some problems when applied to high-dimensional
systems, and present a modified form that solves these problems.

We use the annulus to illustrate the method by examining some of the shadowing
properties of our annulus model using the perfect model scenario (PMS). The
aim is to demonstrate that our method yields shadowing times that are
reasonable. For example, we build a statistical distribution showing how long
a candidate might expect to shadow observations if it starts a particular
distance from the system’s true state. When the initial conditions lie off the
system’s attractor then a reasonable method for identifying shadowing times in
the PMS should, in general, yield longer times the closer the initial
condition is to truth. This result will establish a benchmark for future work
when we search for shadowing trajectories by generating candidates from
observations using the gradient descent method (Judd, 2003; Judd et al.,
2004). As a sanity check we also use a simple distance-based metric to
determine whether our shadowing times are reasonable; they should capture the
time the model visually diverges from the observations. The PMS provides a
robust framework for investigating shadowing in complex systems, as it allows
many aspects of the setup to be controlled, and it also allows us to avoid
dynamical problems due to the $\epsilon$-shadowing properties of the model.

In Sect. 2 we present an analysis of the Smith et al. (2010) shadowing
definition and describe our alternative for high-dimensional systems. Then in
Sect. 3 we describe the rotating annulus model and the experimental setup,
before investigating how the shadowing time depends on the distance of the
candidate state from truth in Sect. 4. Finally we discuss our results and draw
conclusions in Sect. 5.

##  2 A new method for measuring shadowing times in high-dimensional systems

###  2.1 The perfect model scenario

Consider a model $f$ of a system $\tilde{f}$. If the model is a perfect
representation of the system, $f\equiv\tilde{f}$, then there always exists a
trajectory of the model that produces an infinitely long perfect shadow of the
system. To investigate the limits of shadowing it is useful to work in the
perfect model scenario (PMS), where the model and the system are equivalent
($f\equiv\tilde{f}$), but noise does obscure the system state. In the PMS, the
only limit to shadowing the system is the uncertainty in a trajectory’s
initial state. When modelling physical systems this is never the case, because
models are always imperfect representations of the system
($f\not\equiv\tilde{f}$). We call this setting the imperfect model scenario
(IMS), and then model imperfections and dynamical aspects of shadowing a
numerical solution (see below) limit the shadowing capabilities of the model.
In practice, in real systems model imperfections will be the overwhelmingly
larger of the two.

In this paper we only use the PMS; the main advantage being that we know the
true state of the system. This allows us to generate observations by purposely
obscuring the true state by a known noise model of observational error
${\bf\Sigma}$. Controlled experiments can then provide insights into the
properties of shadowing trajectories under initial condition uncertainty
alone. We measure how long a model trajectory started a given distance from
the system state $f({\bf x}+{\bf\delta})$ stays consistent with observations
of the system $\tilde{f}({\bf x})+{\bf e}$, where ${\bf e}\sim{\bf\Sigma}$
represents a realisation of the observational noise model. In the future, our
aim is to explore the IMS as well by allowing the model to be imperfect; Smith
(2002) discusses the subtleties of the PMS and IMS in more detail.

There is a rich literature on the shadowing of numerical solutions of
differential equations by their true solutions, the so-called
$\epsilon$-shadowing problem (Hammel et al., 1987; Grebogi et al., 1990; Sauer
and Yorke, 1991; Sauer et al., 1997; Lai et al., 1999, for example).
Representing differential equations on a computer introduces round-off error
into the numerical solution at each timestep, which tends to grow
exponentially when the system is chaotic until the analytical and numerical
solutions diverge (Hammel et al., 1987). In the current context the PMS
guarantees that this problem does not arise. By construction it is the
numerical solution itself that is the system $\tilde{f}$, not the underlying
equations. Hence the map from time $t$ to $t+\delta t$ is the same for both
model $f$ and system $\tilde{f}$, both including the same numerics for a given
$\bf x$, and so by construction we avoid the $\epsilon$-shadowing problem
entirely. Even when shadowing the laboratory experiment this phenomenon will
not be relevant, because model imperfections will be the dominant source of
error, leading to model drift (see, e.g., Orrell et al., 2001, for an example
in weather modelling).

###  2.2 Consistency in shadowing

The candidate shadowing time $\tau_{S}$ for a specific candidate trajectory is
the time that trajectory remains consistent with observations, given a model
of observational error. Smith et al. (2010) define the two as consistent if
the distribution of the forecast residuals, ${\bf e_{t}\sim{\bf\Phi}}$, and
the distribution of noise realisations given by the noise model,
${\bf\tilde{e}}\sim{\bf\Sigma}$, are equivalent. As long as
${\bf\Phi}\equiv{\bf\Sigma}$, the candidate trajectory is said to shadow the
observations. For a particular model and observations, the shadowing time for
that pair is the maximum among all possible candidates:
$\tau_{S_{t}}=\max_{\bf x}\tau_{S}({\bf x},t)$. By measuring $\tau_{S}$ for a
finite set of candidate trajectories one can only calculate a lower bound on
$\tau_{S_{t}}$.

In practice ${\bf\Phi}$ is expected to differ from ${\bf\Sigma}$ and we aim to
quantify how similar the two distributions are. The degree of similarity
between ${\bf\Phi}$ and ${\bf\Sigma}$ can be assessed by comparing equivalent
percentiles of the two distributions. Suppose we have a single sample
$x_{1},\ldots,x_{N}$ from a distribution ${\bf\Phi}$ and want to determine
whether $\bf\Phi$ is close to some known distribution $\bf\Sigma$. We can
proceed by computing the $r$-th order statistic $X_{(r)}$ of the sample from
$\bf\Phi$ and perform a statistical test to determine whether $X_{(r)}$ is
close to the $r$-th order statistic $\tilde{X}_{(r)}$ from a sample
$\tilde{x}_{1},\ldots,\tilde{x}_{N}$ from $\bf\Sigma$. This calculation
requires the distribution $\tilde{F}_{(r)}$ of the order statistic
$\tilde{X}_{(r)}$. The intuitive way to do this is to repeatedly sample $N$
values $\tilde{x}_{1},\ldots,\tilde{x}_{N}$, from $\bf\Sigma$, calculate
$\tilde{X}_{(r)}$ for each sample, and construct the empirical distribution
$\tilde{F}_{(r)}$. We can then check whether $X_{(r)}$ falls within a
particular confidence interval of $\tilde{X}_{(r)}$. If it does, we say that
the order statistic $X_{(r)}$ is consistent with the distribution
$\tilde{F}_{(r)}$ of the order statistics $\tilde{X}_{(r)}$ and hence the two
distributions are sufficiently close.

Figure 2: Schematic of the general order statistics approach to calculating a
candidate’s shadowing time. The top panel shows a histogram (grey) of the
values generated from a sample of size 1000 from $N(0,1)$. The black
histograms show, for 500 such samples, the distributions of the 50th and 90th
percentiles of each sample. These percentile distributions are shown in the
middle (50th) and lower (90th) panels. On the lower panels vertical lines are
added at percentiles 10% and 90% in those distributions. The 50th and 90th
percentiles of a single test distribution drawn from $N(0,1.01)$ are shown as
grey circles.

For a more robust assessment more than one order statistic $X_{(r)}$ can be
tested in this way, which we illustrate in Fig. 2. In that figure we draw 500
samples of size 1000 from a known noise distribution ${\bf\Sigma}\sim N(0,1)$,
and calculate the 50th ($\tilde{X}_{(50\%)}$) and 90th ($\tilde{X}_{(90\%)}$)
percentiles of each sample. The black bars in the figure show the
distributions of these percentiles. We then produce a test sample of the same
size but from a slightly different distribution ${\bf\Phi}\sim N(0,1.01)$, and
calculate its order statistics $X_{(50\%)}$ and $X_{(90\%)}$, shown by circles
in the lower figures. To retain the null hypothesis that the test sample is
drawn from the same distribution as the noise model at significance level
$p=0.2$, the two circles must fall within the vertical lines enclosing the 80%
confidence intervals of the percentile distributions. In this particular
example they do not, so the two distributions are not sufficiently close.

Applying this procedure to a sample of residuals ${\bf e}_{t}$ and
realizations of a noise model ${\bf\tilde{e}}_{t}$ determines whether the
model trajectory is sufficiently close to shadow the target at a given
significance level.

###  2.3 The Smith et al. (2010) method for measuring shadowing times

Smith et al. (2010) use the the Euclidean length of the residual vector
$\|{\bf e}_{t}\|$ as the test statistic to measure the shadowing time in what
we shall refer to as the “sequence” method. This is the distance in the
$N$-dimensional model state space between the model trajectory and the target.
They calculate this quantity at each time and form a distribution of distances
over a trajectory time series, before calculating the median and 90th
percentile of that distribution.

To obtain a sampling distribution of a given order statistic of the noise
model, they draw a sample of noise realizations ${\bf\tilde{e}}$ of size $N$
from $\bf\Sigma$ at each time. They calculate the Euclidean distance
$\|{\bf\tilde{e}}\|$ for each time, giving a distribution of noise distances,
and then locate the median and 90th percentile of this distribution. Repeating
this sampling procedure $M$ times yields a distribution of $M$ medians and $M$
90th percentiles under the noise model. To determine whether the model
trajectory shadows the target at this time, they test the following null
hypothesis: for both order statistics (the 50th and 90th percentiles) the
percentile of the residual distribution is drawn from the distribution of
samples from the noise distribution. The comparison of medians tests the
residual distribution for bias, and the comparison of the 90th percentiles
tests it for dispersion. They use confidence intervals of 90% for the median
and 99% for the 90th percentile. The procedure begins by using the trajectory
segment with $t=1$ only. If this segment shadows the target then they repeat
with the trajectory segment increased in length by one, and so on, increasing
the trajectory segment by one point at a time. The candidate shadowing time is
determined once a trajectory segment rejects the null hypothesis.

####  2.3.1 The need for an alternative method

We see three problems with the Smith et al. (2010) definition that are
particularly important when $N$ becomes large.

First, Smith et al.’s (2010) sampling approach is computationally expensive.
At each time one requires $M$ samples of size $N$ from the noise model to
produce the distribution of percentiles. Of course $M$ should be as large as
possible to ensure the best correspondence to the true underlying distribution
of each percentile. For systems of very low dimensionality this is feasible
but in higher dimensions this sampling approach is no longer computationally
practical.

Second, a single outlier at one point in the model state space in a very short
trajectory may distort the distribution of the $\|{\bf e_{t}}\|$ and hence
cause the model to fail to shadow for that trajectory length, while a longer
trajectory including that point would reduce the effect of that outlier in the
distribution and hence that longer time series might then shadow.

Finally, the length of the trajectory is important for another reason when
using a distance method. If the trajectory is only length five, for example,
then the distribution of distances is too small to obtain a sensible value of,
say, the 90th percentile of that distribution. For a sensible shadowing time,
therefore, the trajectory must be sufficiently long for statistical accuracy,
but often the model is not capable of shadowing for the minimum of, say,
$O(100)$ leadtimes required.

###  2.4 The “state” method for measuring shadowing times

To avoid some of these problems, we now present a slightly modified version of
the method used to measure the shadowing time from the consistency between
model trajectory and observations that is better suited to higher-dimensional
systems. We also believe this method is more intuitive than using a
distribution of distances from model to target along a trajectory. We call
this the “state” method.

####  2.4.1 Avoiding sampling problems in high-dimensions

First we present a theoretical result that overcomes the sampling
impracticality of the sequence method when applied to cases with $\gg O(10)$
dimensions, opening it up for use in the higher-dimensional context.

Following David (1981), assume independently and identically distributed (IID)
samples of size $N$ from a cumulative distribution function (CDF) $P(x)$, each
sample providing $N$ random variables $X_{1},\ldots,X_{N}$. When the sample
elements are sorted in ascending order, i.e. $X_{(1)}\leq X_{(2)}\ldots\leq
X_{(N)}$, we call $X_{(r)}$ the $r$th order statistic. Denote by $F_{(r)}$ the
CDF of the $r$th order statistic $X_{(r)}$. David (1981, pp. 8–9) showed that
$F_{(r)}$ is

$\displaystyle F_{(r)}(x)$ $\displaystyle=$ $\displaystyle Pr\\{X_{(r)}\leq
x\\}$ (1) $\displaystyle=$ $\displaystyle\sum^{N}_{i=r}{N\choose
i}P^{i}(x)\left[1-P(x)\right]^{N-i}$ (2)

The relation between the binomial sum and the incomplete $\beta$-function
(NIST, 2010, Eq. 8.17.5) means we can replace Eq. (2) with

$F_{(r)}(x)=I_{P(x)}(r,N-r+1)$ (3)

Suppose we are interested in the 90th percentile of the distribution of the
$r$th order statistic $X_{(r)}$ (i.e. $F_{(r)}=0.90$), with $X$ being sampled
from $P(x)$. We require the $x$ that corresponds to the 90th percentile of the
distribution of $X_{(r)}$. Since $F_{(r)}=I_{P(x)}(r,N-r+1)$, we may find the
probability that $P(x)$ assigns the required value $x$ by taking the inverse
of the incomplete beta function, i.e.

$P(x)=I^{-1}_{F_{(r)}}(r,n-r+1)$ (4)

Finally we obtain the value $x$ which gives the required percentile of the
distribution of $X_{(r)}$ using the inverse of $P(x)$: $x=P^{-1}(x)$.

This theoretical result provides an exact distribution of any order statistic
$X_{(r)}$ at negligible computational cost. We have overcome the
intractability of sampling in high dimensions, and so under the assumption of
IID samples we can now perform formal hypothesis testing on a sample of any
size.

####  2.4.2 Definition and properties of the “state” method

With this result in mind, we present a modified form of the Smith et al.
(2010) sequence method used to calculate shadowing times. Our method overcomes
the problems listed in Sect. 2.3.1 above while remaining computationally
practical.

At each time $t$ in a model trajectory we have a model state ${\bf x}_{t}$, a
target state ${\bf\tilde{x}}_{t}$ (whose observational error is well-defined
by the distribution $\bf\Sigma$ in the PMS), and a residual error vector ${\bf
e}_{t}={\bf x}_{t}-{\bf\tilde{x}}_{t}$, each with dimension $N$. Starting at
$t=0$, we test the null hypothesis at significance level $p$ that the residual
error vector ${\bf e}_{t}$ is a sample drawn from the observational error
distribution $\bf\Sigma$. To perform the test we measure the median and the
90th percentile of the residual error vector, and find the distributions of
the equivalent percentiles of the noise distribution analytically using the
order statistics method discussed above. To shadow the observations at time
$t$ we require both the median and 90th percentiles of the residual
distribution to fall within the respective confidence intervals of the noise
distribution at significance level $p$. If the model shadows at $t=0$ then we
repeat the procedure for $t=1,2,\ldots$, and the shadowing time for that model
trajectory is the longest time for which the null hypothesis is retained.

This method differs from the “sequence” method in a number of ways. First, we
avoid sampling the noise distribution by using the inverse $\beta$-function as
described in the previous section222Note that now this result has been
obtained, sampling can be avoided for the “sequence” method as well.. This
avoids problem one in Sect. 2.3.1. Second, our method applies to each point in
the model trajectory in turn, rather than using the whole trajectory at once.
Our method uses the distribution of values over the components of the state
vector, while Smith et al. (2010) map the vector of components to a single
number. This property means we avoid problem two in Sect. 2.3.1. Finally, we
avoid problem three in Sect. 2.3.1 by requiring the state vector dimension $N$
to be large. Therefore our method is not appropriate for use with low-
dimensional systems, where $N$ is not large enough to compute sensible
percentile statistics for the residual error distribution. Take, for instance,
the Lorenz (1963) system, with $N=3$. Determining the 90th percentile in a
sample of three values does not produce a sensible result. Our method is
suitable for medium- and high-dimensional systems, with a dimensionality of at
least $N=O(100)$ and preferably $N=O(\gg 100)$.

We should note that a caveat with our method is that by examining each
trajectory point in turn we may not detect small but systematic outliers. If
an outlier is small but present at each trajectory point the model should
shadow at that point. If the whole trajectory is considered, however, as in
Smith et al. (2010), the outliers at all times will accumulate and the failure
of the model trajectory will be detected.

##  3 Experimental methods

To demonstrate our method for measuring shadowing times from the previous
section and investigate some shadowing properties of our rotating annulus
model we run some perturbation experiments; this section describes how they
are set up.

###  3.1 The rotating annulus simulation

Figure 1 shows a schematic of the rotating annulus with the particular
configuration used in this work, which is the “small annulus” setup used by
Hignett et al. (1985). We work in the PMS, so ignore details of the laboratory
setup that introduce complications due to imperfections in the apparatus, for
example, and we take the numerical representation of Fig. 1 as a perfect
representation of the system.

The mathematical model used to simulate the rotating annulus experiment is the
Met Office / Oxford Rotating Annulus Laboratory Simulation (MORALS) (Farnell
and Plumb, 1976; Hignett et al., 1985; Read et al., 2000). This model solves
the Navier-Stokes, heat transfer, and continuity equations subject to the
Boussinesq approximation, in cylindrical polar coordinates (see Appendix).
There are four prognostic variables: three velocity directions $\bf u$
(radial), $\bf v$ (azimuthal), $\bf w$ (vertical), and temperature $\bf T$.
The fifth field is a kinetic pressure field ${\bm{\Pi}}\equiv{\bf
p}/\rho_{0}$, which is diagnostic (calculated using a Poisson equation from
the other four fields); $\rho_{0}$ is the fluid density at a reference
temperature. The fields are discretized on a staggered Arakawa C grid (Arakawa
and Lamb, 1977), and are non-uniform in the radial and vertical directions to
resolve the boundary layers. All velocities are set to zero at the boundaries,
and the temperature gradient is zero across the top and bottom boundaries. The
parametrisation of fluid properties used to set up a 17% glycerol / 83% water
fluid (by volume) is described in Hignett et al. (1985, Table 1). The model
runs with a timestep of $\delta t=0.02$ s and a spatial resolution of
$(N_{R},N_{\theta},N_{z})=(16,32,16)$ nodes in the radial, azimuthal, and
vertical directions respectively. The number of independent variables333From
$16\times 32\times 16\times 4$ minus the number of points on and outside the
fluid boundary (some outside points are required by the boundary conditions
when using a staggered grid). is $N=24192$.

(a) Snapshot of an equilibrated wavenumber-3 temperature structure in a
horizontal section at $z=5.61$ cm after 3000 s.

(b) Time series of the wavenumber-3 temperature mode amplitude in the same
simulation. The amplitude is found at each timestep from a Fourier transform
of the temperature field at mid-radius / mid-height.

Figure 3: Temperature output from a rotating annulus simulation using the
model described in the text with the setup in Fig. 1.

Fig. 3 shows temperature output from an annulus simulation using this
configuration, with $\Omega=1$ rad s-1 and $\Delta T=4$ degC. Figure 3a is a
typical horizontal temperature field near mid-height, showing a wavenumber-3
structure to the flow, and Fig. 3b shows the amplitude of the wavenumber-3
azimuthal temperature mode as a function of time. The flow starts from rest in
the rotating frame at $t=0$, so there is some transient behaviour in the first
$\sim$1000 s.

To determine whether the flow regime is chaotic we estimated the largest
Lyapunov exponent for this run. Chaos requires aperiodic behaviour in a
deterministic system with at least one positive Lyapunov exponent (Strogatz,
1994, p. 323). To estimate the largest Lyapunov exponent we followed the
procedure that Young and Read (2008) used on annulus simulations in various
flow regimes (Wolf et al., 1985; Hegger et al., 1999), in this case working
with the wavenumber-3 temperature time series. The procedure was the same as
documented in Young and Read (2008), with 62 s the optimal delay time in this
case, and we extended the time series in Fig. 3b to 100000 s to minimise noise
in the calculation. We found the largest Lyapunov exponent to be $0.0061\pm
0.0025{\rm\,bits\,\,s^{-1}}$, which corresponds to an error-doubling time of
1.5–3.0 oscillation periods (125–250 s). This exponent is comparable to values
in the chaotic structural vacillation flow regime
(0.0013–$0.0055{\rm\,bits\,\,s^{-1}}$), and is an order of magnitude larger
than values in non-chaotic flow regimes
(-0.00012–$0.00063{\rm\,bits\,\,s^{-1}}$) (Young and Read, 2008), so we might
expect perturbations to grow at a similar rate to structurally vacillating
flow. The model is deterministic, and the flow is aperiodic for the length of
the time series, so the positive Lyapunov exponent demonstrates (for all
practical purposes) that this flow is chaotic444The reader should note that in
the laboratory the flow under this external forcing is usually regular
wavenumber-3 or 2, which is also the case when running MORALS at higher
resolution (Young and Read, 2008, Fig. 3b). Because we consider only the PMS
in this work, however, direct correspondence with the laboratory experiment is
not so important here..

###  3.2 MORALS perfect model setting and observations

Let ${\bf x}[t]$ and ${\bf\tilde{x}}[t]$ be model and true system states at
time $t$, and $f^{\Delta t}$ and $\tilde{f}^{\Delta t}$ be integration of
states by the model and system for time $\Delta t$. So ${\bf x}[t+\Delta
t]=f^{\Delta t}({\bf x}[t])$ and ${\bf\tilde{x}}[t+\Delta t]=\tilde{f}^{\Delta
t}({\bf\tilde{x}}[t])$ describe the model and system exactly. Let ${\bf
s}[t]={\bf\tilde{x}}[t]+\tilde{\bf e}$ be observations of the system at time
$t$, and $\Delta{\bf x}$ be perturbations to initial states. In the PMS we
simply choose the system being modelled, $\tilde{f}$, to be MORALS with the
same parameters as the model: $f\equiv\tilde{f}$.

We generate a true system trajectory and observations of that trajectory as
follows. First we run MORALS for 1500 s to spin up the model from rest until a
coherent flow structure is present (e.g. Fig. 3) and any transient behaviour
has decayed. From this we run MORALS for an additional 1500 s to obtain an
estimate, at each grid point, of the range of values admitted by the model
when it is in dynamical equilibrium. This gives the natural variability of the
system at each point in space. Let the vector ${\bf r}$ be this range, where
each element contains the range of values that encloses 99% of the values over
these 1500 s at that grid point:

${\bf
r}={\bf\tilde{x}}[t_{0}-1500s,...,t_{0}]_{99.5\%}-{\bf\tilde{x}}[t_{0}-1500s,...,t_{0}]_{0.5\%}$
(5)

The final state in this run (i.e. after a total of 3000 s) will be the first
state in the true system trajectory. To make things simpler later, we redefine
this point as $t=0$. We then obtain a long system trajectory by integrating
this state using MORALS for time $T$ to obtain
${\bf\tilde{x}}[t],\,t\in[0,\ldots,T]$.

We then generate observations ${\bf s}[t]$ every 5 s along the system
trajectory by adding a vector of random numbers $\tilde{\bf e}$ (noise) to the
system trajectory. The random numbers are independently and identically
distributed (IID) and are originally drawn from the Gaussian distribution
$N(0,1)$. These random numbers are converted to observational error by
multiplying by a fixed fraction $\sigma$ ($=0.1$ throughout this paper) of the
natural variability $\bf r$ at each grid point, i.e.

$\tilde{\bf e}\sim\sigma\,{\bf r}\,N(0,1)$ (6)

The effect of the scaling by $\bf r$ ensures that, statistically, the amount
of noise added at each point and to each field is the same fraction of the
natural variability. For this reason, much of the analysis to follow will be
done in terms of quantities expressed as a fraction of the natural
variability, as this non-dimensional metric provides a convenient way to
compare and combine the variables. We shall call this state space the “scaled
space”, as opposed to the (dimensional) “physical space”. For two points in
the physical state space, ${\bf x}_{1}$ and ${\bf x}_{2}$, their separation
$\Delta{\bf x}^{\prime}$ in the scaled space is

$\Delta{\bf x}^{\prime}\equiv({\bf x}_{2}-{\bf x}_{1})^{\prime}=({\bf
x}_{2}-{\bf x}_{1})\circ{\bf r}^{-1}\equiv\Delta{\bf x}\circ{\bf r}^{-1}$ (7)

where $\circ$ denotes the Hadamard (pointwise) product.

###  3.3 Setup of the perturbation experiments

To investigate the shadowing properties of the rotating annulus model MORALS
in the PMS, we calculate trajectories of the model starting from initial
candidate states perturbed from the true state at $t=0$. We use these
candidate trajectories to demonstrate the “state” method developed in the
previous section, by calculating how long each candidate shadows the
observations in Sect. 3.2. By varying the size of the perturbation we shall
determine how the distribution of candidate shadowing times depends on the
original distance of the candidate from the true state.

To generate a candidate state perturbed from the true state
${\bf\tilde{x}}[0]$ we add a random perturbation $\Delta{\bf x}$ as follows.
First, we draw a random vector $\bf n$ of length $N$ from the normal
distribution $N(0,1)$. Second, we scale the Euclidean length of the vector to
be $\delta$, i.e.

${\bf n}^{\prime}=\delta\frac{\bf n}{\|{\bf n}\|}$ (8)

Finally, to obtain a perturbation in physical space (i.e. the perturbation
that is actually added to generate the candidate state), we multiply ${\bf
n}^{\prime}$ by the natural variability at each grid point, $\bf r$. Hence the
perturbation $\Delta{\bf x}$ is

$\Delta{\bf x}=\delta\frac{\bf n}{\|{\bf n}\|}\circ{\bf r}$ (9)

By construction,

$\|\Delta{\bf x}\circ{\bf r}^{-1}\|=\delta$ (10)

so $\delta$ is the length of the perturbation in the scaled space. All the
perturbations of length $\delta$ lie on a $N$-dimensional hypersphere of
radius $\delta$ centred on the truth, and the Gaussian distribution guarantees
that the perturbations are distributed uniformly in direction. By repeating
this procedure any number of candidates may be generated.

This perturbation method is rarely used to efficiently sample initial
condition uncertainty; usually methods such as breeding vectors (Toth and
Kalnay, 1993), singular vectors (Buizza and Palmer, 1995), or even so-called
“nightmare” vectors (Smith et al., 1999) are used, which all aim to give
shorter candidate shadowing times by sampling the fastest-growing directions.
Here we just use a simple method to demonstrate our method for calculating
shadowing times. We note that the distribution of candidate shadowing times
would probably be different were we to use these other methods.

###  3.4 Measuring candidate shadowing times

For each candidate state used to begin a model trajectory we calculate the
candidate shadowing time $\tau_{S}$ with respect to the observations using the
“state” method described in the previous section. At each time $t$ we use the
model $f^{t}({\bf x}[0])$ and observations ${\bf s}[t]$ to construct the
residual error vector ${\bf e}[t]=f^{t}({\bf x}[0])-{\bf s}[t]$. The
observational error $\tilde{\bf e}$ is well-defined by Eq. (6). We determine
whether the distribution of scaled forecast residuals ${\bf e}[t]\circ{\bf
r}^{-1}$ is consistent with the distribution of scaled observational error
$\tilde{\bf e}\circ{\bf r}^{-1}$, testing the null hypothesis at significance
level $p$ for the 50th and 90th percentiles. If either test statistic lies
outside the confidence interval of the noise distribution then the candidate
trajectory does not shadow the observations at time $t$. The candidate
shadowing time is the final $t$ at which the test is passed, counting upwards
from zero.

The shadowing time should capture the time the model diverges from the
observations. Hence as a sanity check on our shadowing times we also use a
simple measure of how long the model remains close to observations. The
Euclidean norm $\|{\bf x}\|=({\bf x}^{\prime}{\bf x})^{1/2}$ of a vector of
length $N$ whose components are drawn from the Gaussian distribution
$N(0,\sigma^{2})$ follows a $\chi$-distribution with $N$ degrees of freedom
and scale parameter $\sigma$ (Cramér, 1946, p. 236). For $N\to\infty$ the
asymptotic form of the $\chi$-distribution is a normal distribution with mean
$\sigma\sqrt{N-1/2}$ and variance $\sigma^{2}/2$ (Stuart and Ord, 1994, p.
544–5). So if a vector of $N$ observations $s_{i}$, $i=1,\ldots,N$ is
constructed by adding IID random numbers with distribution $N(0,\sigma^{2})$
to a true state $x_{i}$, $i=1,\ldots,N$, then the expected Euclidean distance
from the truth to the observations in our scaled space is very close to
$\sigma\sqrt{N}$ for the value of $N$ we use.

The Euclidean distance in the scaled space between a model candidate
trajectory $f^{t}({\bf x}[0]))$ and the observations ${\bf s}(t)$ is

$D(t)=\|\\{f^{t}({\bf x}[0])-{\bf s}[t]\\}\circ{\bf r}^{-1}\|$ (11)

so we say that the candidate trajectory remains close to the observations as
long as

$D(t)\leq m\sigma\sqrt{N}$ (12)

is satisfied, i.e. once the distance from model to observations becomes
greater than $m\sigma\sqrt{N}$ then the candidate trajectory is no longer
close. We call this time $\tau_{D}$. Based on test data we use $m=2$, as this
represents a considerable deviation from the expected distance between
observations and truth, but it turns out that $\tau_{D}$ is not very sensitive
to $m$. Now, this is basically a simple RMS error metric, and while that is
not an optimal error measure for nonlinear models (see, e.g., McSharry and
Smith, 1999), the purpose of using this metric is to catch when the trajectory
is not close to the observations, and the RMS error will be large when that is
the case, so it is useful as a secondary measure here.

##  4 Results

We are interested in how long a candidate trajectory might expect to shadow
observations if the setup of an experiment is perfect and it starts at a given
distance from the truth. This is an important benchmark to measure, because in
practice with physical systems the setup is never perfect. A shadowing time
measured in the imperfect model scenario or using laboratory data is of
limited use if we have nothing to compare it against. By using the PMS we can
obtain an upper limit to the shadowing properties of a particular system under
this type of initial condition uncertainty. We expect a priori that, for
perturbations in random directions, the time a trajectory remains close to
observations should decrease as its initial distance from truth increases.

###  4.1 Candidate trajectories started from perfect initial conditions

First, however, we must calibrate the “state” method by selecting an
appropriate significance level $p$. To do so we ran a “null” candidate
trajectory started from perfect initial conditions, i.e. with $\delta=0$. We
measured the shadowing time against 32 different realisations of the
observational noise model. By construction, the real shadowing time for this
case is infinity, but an intrinsic property of our “state” method (and any
other test based on repeated significance tests, including the Smith et al.
(2010) “sequence” method, is that some fraction of the candidates will fail to
shadow by a given time, even if the forecast is perfect, because the method is
based on significance testing. This fraction depends critically on the
significance level at which the null hypothesis is tested. So by calculating
the parameters required to avoid this problem we can determine a priori the
minimum significance level at which our “state” method should be used in
practice.

Figure 4: Consequences of the Type I error property of the “state” method for
low significance levels. For each significance level (indicated by its
$p$-value) the solid line shows the expected fraction of candidates (from Eq.
13) that fail to shadow because of a Type I error as a function of the number
of verifications, and the dotted line shows the actual fraction that failed
out of the 32 noise realisations. The dotted lines stop at 300 verifications,
corresponding to maximum trajectory lengths of 1500 s.

The problem occurs because at each verification time along the trajectory we
compute two test statistics, both of which must fall within prescribed
confidence intervals in order to retain the null hypothesis at significance
level $p$ and hence shadow at that time. The probability of falsely rejecting
the null hypothesis when it is true (Type I error) is therefore $1-(1-p)^{2}$.
Within the first $n$ verifications (in this case we check the shadow every 5 s
along the trajectory), the probability of at least one Type I error is

$r=1-(1-p)^{2n}$ (13)

which is the expected fraction of candidates that will fail to shadow within
the first $n$ verifications purely from Type I errors. For significance levels
commonly used in statistical testing ($p=0.1$, 0.01, say) this probability
becomes large after $O(10)$ verifications. We must avoid these Type I errors
because they erroneously assign a candidate shadowing time shorter than the
true value.

In Fig. 4 we compare the fraction of candidates that failed to shadow because
of a Type I error at each $t$ against theoretical predictions from Eq. (13);
the two correspond well. By $t=50$ s ($n=10$), at significance level $p=0.1$
almost all the candidates have rejected the null hypothesis and hence recorded
a (false) shadowing time, as have almost all of the candidates by $t=500$ s
($n=100$) for $p=0.01$. The 25th percentile of the distribution for $p=0.001$
also falls around this time, there is a sole outlier for $p=0.0001$, and for
$p\leq 10^{-5}$ none of the candidates suffer a Type I error.

This property is a fundamental limitation of any method that relies on
repeated statistical significance tests, but we can largely avoid it in
practice through our choice of $p$. We can find the largest value of $p$ such
that in the event that the true shadowing time is the equivalent of $n$
verifications, fewer than $R$ out of a total of $E$ candidates will suffer a
Type I error. This condition is $r<R/E$, and hence

$p<1-\left(1-\frac{R}{E}\right)^{1/2n}$ (14)

is required, by rearranging Eq. (13). This problem with Type I errors is well-
known in the statistical and medical literature (Abdi, 2007); the expression
in Eq. (14) is an example of a Šidák equation.

###  4.2 Candidate trajectories started a fixed distance from truth

Now we demonstrate our method for identifying shadowing times by examining how
long candidate trajectories started a fixed distance from truth shadow the
observations.

(a) Temperature model trajectory (black line) and corresponding observations
(grey line).

(b) Temperature residual error (black line) and corresponding $\pm\sigma$ and
$\pm 2\sigma$ observational error (dotted lines).

Figure 5: Example annulus candidate trajectory showing the temperature time
series and residual error at a single grid point: $R=4.81$ cm, $\theta=6.19$
rad, $z=5.61$ cm. This particular candidate trajectory begins a distance
$\delta=0.1$ from the truth in the scaled space (Eq. 10).

We generated eight sets of initial candidate states using the method in Sect.
3.3, seven sets containing 64 candidates and one containing 256 (in order to
have one with better statistics). While each set may be small compared with
$N$, because the perturbation directions are random and initial conditions
almost surely lie off the attractor, we would simply expect a smoothing of the
distribution of candidate shadowing times if the number of candidates were
increased further. Much longer or shorter shadowing times might occur in the
special case of a perturbed initial condition falling on the attractor, but
this occurs with probability zero. Within each set the distance from truth of
each candidate was the same in the scaled space. We used $\delta=10^{-4}$,
$10^{-3}$, $10^{-2}$, $e.10^{-2}$, $e.10^{-1}$, $1$, and $e$ as the initial
distances from truth in the scaled space for the 64-candidate sets, and
$\delta=0.1$ for the set with 256 candidates555Multiples of $e$ were used to
ensure approximately regular spacing of the perturbation sizes on a
logarithmic scale.. Using $\sigma=0.1$ the expected distance between
observations and truth is $0.1\sqrt{24192}=15.6$ (Sect. 3.4). Hence the
largest distance between candidate state and truth ($\delta=e$) corresponds to
about 15–20% of the observational error, $\delta=0.1$ corresponds to about
0.6%, and $\delta=10^{-4}$ corresponds to about 0.0005% of the observational
error. The candidate trajectories were integrated by MORALS for 1500 s of
simulated time, recording the state every 5 s.

To calculate the appropriate $p$-values to use when measuring those
candidates’ shadowing times we used Eq. (14) with either $E=64$ or $E=256$ and
$n=301$ (maximum trajectory length of 1500 s, verifying every 5 s, plus a
verification at $t=0$). So for $R<1$ (fewer than one candidate suffering a
Type I error) we found that $p<2.6\times 10^{-5}$ when $E=64$ and $p<6.5\times
10^{-6}$ when $E=256$, so we used $p=10^{-5}$ for the 64-candidate sets and
$p=10^{-6}$ for the 256-candidate set.

Figure 5a shows one grid point in a typical MORALS candidate trajectory and
corresponding (temperature) observations in the PMS using the setup described
in Sect. 3.1. The residual error between the trajectory and observations is
shown in Fig. 5b, and the trajectory in Fig. 5a diverges from the observations
at this particular grid point around 800 s, corresponding to the time the
residual in Fig. 5b grows, on average, to be larger than the $\pm 2\sigma$
error bars.

We measured the distance between the trajectories and observations as a
function of time:

$D(t)=\|\\{f^{t}({\bf x}[0])-{\bf s}[t]\\}\circ{\bf r}^{-1}\|$ (15)

Figure 6 shows this distance for each trajectory in the 256-candidate set (an
initial distance from truth of size $\delta=0.1$ in the scaled space). By
construction, in the scaled space the expected initial distance between each
candidate and the observations is $\sigma\sqrt{N}\sim 15.6$.

Figure 6: Distance between 256 candidate trajectories and one realisation of
the observational noise (using Eq. 15) for candidates begun a distance
$\delta=0.1$ from truth in the scaled state space introduced in Eq. (7). The
vertical scale shows the Euclidean distance between candidate trajectory and
observations in the scaled state space, and is logarithmic. The dotted line is
at $D=\sigma\sqrt{N}$. This distribution of error growth over the candidates
is robust with respect to different realisations of the observational noise
(four were tried in total, with similar results).

In this figure we see a rapid growth of the distance between candidate
trajectories and observations from the expected initial observational error to
the size of the attractor over a period of $\sim$300 s between $t=500$ s and
800 s. Intuitively, we require our method to show that a candidate shadows
observations until the distance between candidate and observations diverges
from from the expected initial distance $\sigma\sqrt{N}$.

Using the “state” method described in Sect. 3.4, we measured the candidate
shadowing times for each candidate in the eight sets, using a range of eleven
significance levels in each case from $p=0.1$ up to the maximum significance
permitted by floating point precision ($p=10\epsilon\approx 4.4\times
10^{-15}$). Figure 7 shows the distribution of candidate shadowing times for
the trajectories in the 256-candidate set ($\delta=0.1$) for each of these
significance levels.

Figure 7: Distributions of candidate shadowing times measured using the
“state” method for the 256 candidate trajectories in Fig. 6 with $\delta=0.1$.
Distributions for eleven different significance levels ($p$-values) are shown,
with the $p$-value on the $x$-axis. Each box-whisker plot shows the
distribution of candidate shadowing times at one significance level. The box
represents the interquartile range, the middle line is the median, and the two
whiskers represent the maximum/minimum or the 75%/25% values +/- 1.5 times the
interquartile range, whichever is closest to the 75%/25% value. Outliers are
shown with circles.

Although above we calculated the “appropriate” significance level to use in
this context, it is important to determine how sensitively the results depend
on this quantity. The figure shows that our candidate shadowing times are
generally not very sensitive to the significance level. Once the candidate
trajectory diverges from the observations, it does so rapidly, reflecting the
trajectory’s rapid divergence from the observations as shown in Fig. 6. As the
residual error distribution diverges from the observational error
distribution, the test rejects the null hypothesis even at the highest
significance levels, and the robustness of our method in this respect is very
encouraging. Below $p=0.001$ the candidate shadowing times are significantly
smaller, which is expected from the analysis of the “null” trajectories above.

At $p=10^{-6}$, the “appropriate” significance level for this candidate set,
the candidate shadowing times are in the range 450–700 s. This concurs with a
visual analysis of the trajectories in Fig. 5 and the visual estimate of about
500 s from Fig. 6.

To determine whether the shadowing times are sensitive to position on the
model attractor, we compared the candidate shadowing time distributions
against two candidate trajectory sets started from other points in the state
space. To start candidate trajectories from other points in the state space,
the spinup time was changed from 1500 s to 1525 s and 1550 s. At this point in
$(\Omega,\Delta T)$ parameter space the timescale of the fluid’s major
oscillation (which determines the extent of the model attractor) is about 75 s
(Fig. 5a), so trajectories started from these other two spinup times are well-
spaced about the model attractor. At each position we ran two 64-candidate
sets with $\delta=10^{-3}$ and $e.10^{-1}$. The distributions of candidate
shadowing times for $p=10^{-6}$ were very similar, with extrema and medians
varying by only 5-10%. Hence the distribution of candidate shadowing times is
not generally sensitive to position in state space, at least for the points
tested.

Now we have established that our method measures reasonable candidate
shadowing times, we consider what the expected candidate shadowing time should
be given its initial distance from truth. Conversely, how close to the truth
must a candidate trajectory begin at to be expected to shadow for a particular
time? With our perturbation experiments we can approach this question in a
controlled fashion.

Figure 8: Distributions of candidate shadowing times for eight initial
distances from truth $\delta$, using significance level $p=10^{-5}$ for all
but $\delta=0.1$, which uses $p=10^{-6}$. The boxes are constructed in the
same way as in Fig. 7.

Using $p=10^{-5}$ for the 64-candidate sets, and $p=10^{-6}$ for the
256-candidate set, we plot in Fig. 8 the distribution of candidate shadowing
times for the eight different initial distances from truth $\delta$. There is
clearly a logarithmic relationship between these candidate shadowing times and
the initial distance from truth. An approximate empirical best fit line to the
data is

$\tau_{S}\approx 250\left(1-\log_{10}\delta\right)$ (16)

where $\tau_{S}$ is the candidate shadowing time, measured in seconds, and
$\delta$ is the initial distance from truth (Eq. 10). Conversely, we predict
that a candidate will have an expected shadowing time of $\tau_{S}$ when the
initial distance from truth in the scaled space is approximately

$\delta\approx 10^{1-\tau_{S}/250}$ (17)

To test this prediction, we ran a final set of candidate trajectories which we
expect to shadow for around $\tau_{S}=2000$ s. From Eq. (17) this corresponds
to $\delta\approx 10^{-7}$. We started the final 64-candidate set from the
same point in state space as the previous eight sets, with an initial distance
from truth of $\delta=10^{-7}$, and found the distribution of candidate
shadowing times to be close to the predicted 2000 s:

Percentile (%) | 0 | 25 | 50 | 75 | 100
---|---|---|---|---|---
$\tau_{S}$ (s) | 1925 | 2000 | 2090 | 2130 | 2405

Note that the maximum shadowing time $\tau_{S_{t}}$ over all possible
candidates is not what is being measured here. Candidates beginning on the
attractor are special cases that the dynamics allows for but will almost never
occur with our perturbation method. We loosely use the term “on the attractor”
here to include points on the attractor and on the relevant stable manifold.
In these cases a candidate further away from the truth may have a longer
shadowing time, but our perturbation method generates an initial condition on
the attractor with probability zero. Hence we can only conclude that the
expected candidate shadowing time for a trajectory beginning a given distance
from truth varies logarithmically with that distance; the maximum shadowing
time may be much longer.

Our aim was to demonstrate that the method behaves sensibly, and it does so:
the statistical distribution of shadowing times decreases as initial distance
from the truth increases. The functional form of this relationship is close to
logarithmic, however, which is perhaps a little surprising. Smith et al.
(1999) showed that, even for low-dimensional systems like the Lorenz
equations, the error $q$-pling time $\tau_{q}$ behaves such that the error-
quadrupling time $\tau_{4}$ is generally much less than twice the error-
doubling time $\tau_{2}$, indicating explosive and episodic uncertainty growth
(their Fig. 8). By contrast, exponential divergence of two trajectories
implies $\tau_{q^{2}}\approx 2\tau_{q}$. So one would not necessarily expect
an logarithmic dependence of the shadowing time on initial distance from
truth, as that would be equivalent to $\tau_{q^{2}}\approx 2\tau_{q}$. This is
what we have observed, however, in Fig. 8. Our initial perturbations are
possibly too large for this phenomenon to be relevant; Smith et al. (1999) use
infinitesimal perturbations in their analysis. Furthermore, they state that
$\tau_{q^{2}}\approx 2\tau_{q}$ may hold once $q$ is large enough. $q=32$ is
too small for the systems they use, but our experiments represent $q$ as high
as $\approx 400$, which may be sufficient.

###  4.3 Sanity check using the simple distance-based metric

In Sect. 3.4 we introduced a simple distance-based metric to determine the
time $\tau_{D}$ that a candidate trajectory remains within $m\sigma\sqrt{N}$
of the observations in the scaled state space (Eqs 11 and 12). This estimates
the time at which the candidate trajectory diverges from the observations.

Figure 9: Contour plot of the unsigned distance $|{\bf e}[t]\circ{\bf
r}^{-1}|$ between candidate trajectory and observations for each grid point in
a candidate with $\delta=0.1$ (same run as Fig. 5), for each point in the
trajectory. A time series for a single grid point is a single vertical line in
this diagram, and for a particular $t$ all the model grid points are in a
horizontal line (the order of these grid points along the $x$-axis is
unimportant). The horizontal solid line is at $t=\tau_{D}$ for this candidate
(using $m=2$), and the dotted lines show $\tau_{S}$ for the 11 different
$p$-values in Fig. 7. Note that to enhance the contrast, the contour intervals
are irregularly spaced.

The trajectory and residual in Fig. 5 show an example at a single model grid
point, where the model trajectory diverges from the observations around
$t=800$ s, but the definition of $\tau_{D}$ should identify the time at which
this occurs for the whole state vector. In Fig. 9 we show a visualisation of
the same candidate trajectory (of which Fig. 5 was a single grid point) that
extends the visualisation of residual error at one grid point to the whole
state vector. There is an abrupt transition from small distances (dark
regions) to large distances (light regions) over the whole state vector, quite
similar over most of the grid points, which is the time at which the candidate
trajectory diverges from the observations. Below we list several percentiles
of the probability distribution function for $\tau_{D}$ measured for each
trajectory in the 256-candidate set.

Percentile (%) | 0 | 25 | 50 | 75 | 100
---|---|---|---|---|---
Shadowing time (s) | 675 | 735 | 760 | 775 | 980

In Fig. 9, $\tau_{D}$ is identified by a solid horizontal line. There is
general agreement between $\tau_{D}$ and the time at which the absolute
distance between candidate trajectory and observations increases rapidly, i.e.
from dark to light in the figure as we move upwards. $\tau_{D}$ corresponds
well to this transition for all candidates for which such diagrams were
produced, being within $\pm 50$ s in each case. Hence $\tau_{D}$ estimates
well the time at which a candidate trajectory diverges from observations.
Furthermore, the time between diverging from the expected distance between
observations and truth and reaching the saturation distance (distance between
two points randomly chosen on the attractor) is relatively short compared with
the length of a trajectory (Fig. 6). So as long as $m$ is not so large as to
make $m\sigma\sqrt{N}$ larger than the saturation distance, then $\tau_{D}$
will be quite insensitive to $m$.

The dotted lines show $\tau_{S}$ measured in the previous section for this
particular candidate for each of the eleven $p$-values. In this case, as in
the others, $\tau_{S}<\tau_{D}$, which demonstrates that the distribution of
residual error loses consistency with the observational noise distribution
some time before the trajectory and observations diverge visually. Overall,
however, our simple distance metric confirms that our method for measuring
shadowing times is a reasonable measurement of the time model trajectories
diverge from observations.

##  5 Discussion and conclusions

We have presented a study of shadowing in the rotating annulus using the
perfect model scenario. First we developed a version of the Smith et al.
(2010) “sequence” method of calculating the shadowing time that is more
suitable for use with high-dimensional systems. In the course of this analysis
we found a way to avoid the practical sampling problems of applying the
sequence method in higher dimensions. Our “state” method measures the
consistency of the residual error distribution and the observational error
distribution, and requires consistency at each time up to $t$ in order to
shadow for a time $t$.

We demonstrated our state method using rotating annulus simulations in a
series of perturbation experiments. Generating sets of candidate states a
fixed distance from the true state, we measured the candidate shadowing times
with respect to a set of artificial observations. Because our perturbation
method ensured that, with probability one, initial conditions would not lie on
the model attractor or stable manifold, we expected (almost surely) the
trajectories to diverge from the observations, and they did so. The shadowing
times decreased as the trajectories’ initial distances from truth increased,
with a clear logarithmic relationship between the initial distance from truth
and the distribution of candidate shadowing times. We used the relationship to
predict how close to the truth the candidate would have to be for its expected
shadowing time to be 2000 s, and the results agreed well with the prediction.
The candidate shadowing times (Fig. 7) also corresponded well to the time
trajectories remained close to observations measured using a simple distance
metric (Fig. 9), concurred with a visual analysis of the forecasts (Fig. 5),
and appeared to be similarly distributed at other points in the model state
space, all of which strengthen our conclusion that the method we have
presented to measure shadowing times is reasonable.

By using the PMS we were able to measure shadowing times against the true
state of the system using a simple setup. In particular, deliberately
simplifying the problem in this way avoided complications introduced by the
$\epsilon$-shadowing properties of the model, which are currently unknown. We
expect the model’s equations to be non-hyperbolic and so long
$\epsilon$-shadowing trajectories probably don’t exist, for various reasons
(Lai et al., 1999). The larger the number of dimensions the more complications
this is expected to introduce: for example, see Hayes (2003) for an analysis
in $O(100)$ dimensions, which is still a factor of $O(100)$ smaller than
MORALS. To investigate these properties further, one could determine whether
phenomena such as glitching (Sauer and Yorke, 1991; Dawson et al., 1994) or
unstable-dimension variability (Kostelich et al., 1997) might be important.
Dawson et al. (1994) showed that if a finite-time local Lyapunov exponent
fluctuates about zero then $\epsilon$-shadowing will fail quickly, and so
quantifying this using the hyperbolicity exponent (Sauer, 2002), which is
possible in principle although difficult in practice, would provide some
insight. Another simple experiment would be to compare a “clean” model
trajectory with one where noise on the order of the roundoff error is
introduced at each timestep, which would give a lower bound on the
$\epsilon$-shadowing times. This approach has precedent in Sauer et al.
(1997); see also Smith and Spiegel (1987) for an example using numerical
representations of 1–1 chaotic maps without round-off error. In any case, we
expect the model’s intrinsic structural error to be larger than these other
sources of error, so in practice shadowability of the real experiment will
almost certainly be limited by model error rather than anything else.

These results define a benchmark for further study of shadowing in this
system. From Eq. (17) we have a quantitative measure of whether a candidate
shadowing time of, say, 400 s is “good” or “bad”, because it tells us how
close to the true state we should expect a randomly orientated perturbation to
be in order to expect that candidate shadowing time. In addition, in the PMS
the candidate shadowing time tells us approximately how close to the truth our
candidate was initially, assuming random perturbations, and hence how much we
might expect to be able to improve it. Because of our perturbation method our
candidate shadowing times will be lower bounds in the PMS — a method that
generates candidates close to or on the attractor (such as gradient descent)
should yield longer candidate shadowing times. With an imperfect model or
observational data, however, we would expect these shadowing times to be an
upper bound given that then the model and system attractors will then be
disjoint or almost disjoint.

It should be stressed that these experiments were performed with one
particular set of model parameters. The particular candidate shadowing times
obtained would therefore be expected to change for different points in the
parameter space. In a non-chaotic flow regime, for example, it is conceivable
that a candidate trajectory might shadow observations for all $t$.
Nevertheless, in other chaotic regions of parameter space we predict that the
candidate shadowing times should vary in a similar way as a function of
initial distance from the truth.

The work we have presented in this paper falls into a larger body of work
examining the shadowing properties of the rotating annulus system by applying
the method of gradient descent (Judd, 2003). This method searches for
shadowing trajectories by relaxing a sequence of observational states onto the
model attractor. Gradient descent has been used in low- (Judd, 2003, and
elsewhere), intermediate- (Judd et al., 2004) and high- (Judd et al., 2008)
dimensional models for the purposes of state estimation. We intend to combine
our “state” method for measuring shadowing times with an algorithm for
gradient descent and use this to determine whether gradient descent can
produce candidates with long shadowing times in the rotating annulus system.
In the perfect model scenario we will know how far away a state obtained via
gradient descent is from the truth, and hence from Eq. (16) we will have an a
priori estimate of how long we might expect a candidate to shadow
observations. This work will be a useful step towards determining whether the
gradient descent method can produce model trajectories with long shadowing
times in high-dimensional models of weather and climate.

## Acknowledgements

We thank Leonard Smith for providing the original impetus for this work and
for numerous insights during its duration. We thank him and Hailang Du for
reading the manuscript fully more than once. We also thank Kevin Judd, Peter
Read, and Thomas Stemler for useful conversations on a number of topics. We
thank two anonymous reviewers whose comments have greatly helped us to improve
the manuscript. The program used to compute Lyapunov exponents was written by
Peter Read, and the TISEAN package is available from http://www.mpipks-
dresden.mpg.de/$\sim$tisean/. RMBY and FN acknowledge financial support from
the Grantham Research Institute on Climate Change and the Environment. RMBY
also acknowledges financial support from NERC Studentship NER/S/A/2005/13667.

##  Appendix A MORALS equations

MORALS solves the fluid equations in cylindrical polars $(R,\phi,z)$ for a 3D
velocity field ${\bf u}=(u,v,w)$ and temperature $T$ in a reference frame
rotating at angular velocity $\Omega$ using a finite-difference integration
scheme with leapfrog time-stepping. The pressure $\Pi$ is found diagnostically
using a Poisson equation (Farnell and Plumb, 1976, Sect. 5.1). $T$ is relative
to a reference temperature $T_{\rm R}=22$ ∘C, and $\Pi$ is relative to a
reference pressure $\Pi_{0}(R,z)=\frac{1}{2}\Omega^{2}R^{2}+g(d-z)$. The
equations solved are the continuity, heat, and Navier-Stokes equations:

$\nabla\cdot{\bf u}=0$ (18) $\frac{\partial T}{\partial t}+{\bf u}\cdot\nabla
T=\nabla\cdot(\kappa\nabla T)$ (19) $\frac{\partial{\bf u}}{\partial t}+{\bf
u}\cdot\nabla{\bf u}+2{\bf\Omega}\times{\bf u}-(\Omega^{2}R\hat{\bf
R}-g\hat{\bf z})\frac{\Delta\rho}{\rho_{0}}+\nabla\Pi={\bf F}$ (20)

Here $\bf F$ is the viscous term, but is more complicated than the standard
$\nu\nabla^{2}{\bf u}$ because it allows spatial variations of the viscosity
as well; see Farnell and Plumb (1976, Eqs 2.1–2.3, 2.6) for the full
expression. These equations are closed by an equation of state for density
$\rho(T)=\rho_{0}(1+\rho_{1}T+\rho_{2}T^{2})$ with $\Delta\rho=\rho-\rho_{0}$
and two “constitutive” relations for viscosity
$\nu(T)=\nu_{0}(1+\nu_{1}T+\nu_{2}T^{2})$ and thermal diffusivity
$\kappa(T)=\kappa_{0}(1+\kappa_{1}T)$. On all boundaries ${\bf u}={\bf 0}$ and
$\hat{\bf n}\cdot\nabla\Pi=0$, on the top and bottom boundaries $\hat{\bf
n}\cdot\nabla T=0$, and on the inner and outer cylinders $T=T_{a}-T_{R}$ and
$T=T_{b}-T_{R}$ respectively.

## References

  * Abdi (2007) Abdi, H. (2007), The Bonferonni and Šidák corrections for multiple comparisons, in Encyclopedia of Measurement and Statistics, edited by N. J. Salkind, pp. 103–107, Sage.
  * Arakawa and Lamb (1977) Arakawa, A., and V. R. Lamb (1977), Computational Design of the Basic Dynamical Processes of the UCLA General Circulation Model, Meth. Comput. Phys., 17:173–265, doi:10.1016/B978-0-12-460817-7.50009-4.
  * Bowen (1975) Bowen, R. (1975), omega-Limit Sets for Axiom A Diffeomorphisms, J. Differ. Equations, 18:333–339, doi:10.1016/0022-0396(75)90065-0.
  * Buizza and Palmer (1995) Buizza, R., and T. N. Palmer (1995), The singular-vector structure of the atmospheric global circulation, J. Atmos. Sci., 52:1434–1456, doi:10.1175/1520-0469(1995)052<1434:TSVSOT>2.0.CO;2.
  * Cramér (1946) Cramér, H. (1946), Mathematical Methods of Statistics, 11th print ed., Princeton University Press.
  * David (1981) David, H. (1981), Order statistics, second ed., Wiley.
  * Dawson et al. (1994) Dawson, S., C. Grebogi, T. Sauer, and J. A. Yorke (1994), Obstructions to shadowing when a Lyapunov exponent fluctuates about zero, Phys. Rev. Lett., 73:1927–1930, doi:10.1103/PhysRevLett.73.1927.
  * Farnell and Plumb (1976) Farnell, L., and R. Plumb (1976), Numerical integration of flow in a rotating annulus II: three dimensional model, Tech. rep., Occasional Note Met O 21 76/1, Geophysical Fluid Dynamics Laboratory, Meteorological Office, Bracknell, Berkshire.
  * Gilmour (1998) Gilmour, I. (1998), ‘$\iota$-shadowing, probabilistic prediction and weather forecasting’, Ph.D. thesis, University of Oxford.
  * Grebogi et al. (1990) Grebogi, C., S. M. Hammel, J. A. Yorke, and T. Sauer (1990), Shadowing of physical trajectories in chaotic dynamics: Containment and refinement, Phys. Rev. Lett., 65:1527–1530, doi:10.1103/PhysRevLett.65.1527.
  * Hammel et al. (1987) Hammel, S. M., J. A. Yorke, and C. Grebogi (1987), Do numerical orbits of chaotic dynamical processes represent true orbits?, J. Complexity, 3:136–145, doi:10.1016/0885-064X(87)90024-0.
  * Hammel et al. (1988) Hammel, S. M., J. A. Yorke, and C. Grebogi (1988), Numerical orbits of chaotic processes represent true orbits, B. Am. Math. Soc., 19:465–469, doi:10.1090/S0273-0979-1988-15701-1.
  * Hayes (2003) Hayes, W. B. (2003), Shadowing-based reliability decay in softened n-body simulations, Astrophys. J., 587:L59–L62, doi:10.1086/375263.
  * Hegger et al. (1999) Hegger, R., H. Kantz, and T. Schreiber (1999), Practical implementation of nonlinear time series methods: The TISEAN package., Chaos, 9:413–435, doi:10.1063/1.166424.
  * Hide (1953) Hide, R. (1953), Some experiments on thermal convection in a rotating liquid, Q. J. Roy. Meteor. Soc., 79:161, doi:10.1002/qj.49707933916.
  * Hide and Mason (1975) Hide, R., and P. Mason (1975), Sloping convection in a rotating fluid, Adv. Phys., 24:47–100, doi:10.1080/00018737500101371.
  * Hignett et al. (1985) Hignett, P., A. A. White, R. D. Carter, W. D. N. Jackson, and R. M. Small (1985), A comparison of laboratory measurements and numerical simulations of baroclinic wave flows in a rotating cylindrical annulus, Q. J. Roy. Meteor. Soc., 111:131–154, doi:10.1002/qj.49711146705.
  * Ikeda (1979) Ikeda, K. (1979), Multiple-valued stationary state and its instability of the transmitted light by a ring cavity system, Opt. Commun., 30:257–261, doi:10.1016/0030-4018(79)90090-7.
  * Judd (2003) Judd, K. (2003), Nonlinear state estimation, indistinguishable states, and the extended Kalman filter, Physica D, 183:273–281, doi:10.1016/S0167-2789(03)00180-5.
  * Judd et al. (2004) Judd, K., L. Smith, and A. Weisheimer (2004), Gradient free descent: Shadowing, and state estimation using limited derivative information, Physica D, 190:153–166, doi:10.1016/j.physd.2003.10.011.
  * Judd et al. (2008) Judd, K., C. A. Reynolds, T. E. Rosmond, and L. A. Smith (2008), The Geometry of Model Error, J. Atmos. Sci., 65:1749–1772, doi:10.1175/2007JAS2327.1.
  * Kostelich et al. (1997) Kostelich, E. J., I. Kan, C. Grebogi, E. Ott, and J. A. Yorke (1997), Unstable dimension variability: A source of nonhyperbolicity in chaotic systems, Physica D, 109:81–90, doi:10.1016/S0167-2789(97)00161-9.
  * Lai et al. (1999) Lai, Y.-C., C. Grebogi, and J. Kurths (1999), Modeling of deterministic chaotic systems, Phys. Rev. E, 59:2907–2910, doi:10.1103/PhysRevE.59.2907.
  * Lorenz (1963) Lorenz, E. N. (1963), Deterministic nonperiodic flow, J. Atmos. Sci., 20:130–141, doi:10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2.
  * Lorenz (1967) Lorenz, E. N. (1967), Laboratory Models of the Atmosphere, in The Nature and Theory of the General Circulation of the Atmosphere, pp. 114–126, World Meteorological Organization.
  * McSharry and Smith (1999) McSharry, P. E., and L. A. Smith (1999), Better nonlinear models from noisy data: attractors with maximum likelihood, Phys. Rev. Lett., 83:4285–4288, doi:10.1103/PhysRevLett.83.4285.
  * NIST (2010) NIST (2010), NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/.
  * Orrell et al. (2001) Orrell, D., L. Smith, J. Barkmeijer, and T. N. Palmer (2001), Model error in weather forecasting, Nonlinear Proc. Geoph., 8:357–371, doi:10.5194/npg-8-357-2001.
  * Read et al. (1992) Read, P. L., M. J. Bell, D. W. Johnson, and R. M. Small (1992), Quasi-periodic and chaotic flow regimes in a thermally driven, rotating fluid annulus, J. Fluid Mech., 238:599–632, doi:10.1017/S0022112092001836.
  * Read et al. (2000) Read, P. L., N. P. J. Thomas, and S. H. Risch (2000), An evaluation of Eulerian and semi-Lagrangian advection schemes in simulations of rotating, stratified flows in the laboratory. Part I: Axisymmetric flow., Mon. Weather Rev., 128:2835–2852, doi:10.1175/1520-0493(2000)128<2835:AEOEAS>2.0.CO;2.
  * Sauer and Yorke (1991) Sauer, T., and J. A. Yorke (1991), Rigorous verification of trajectories for the computer simulation of dynamical systems, Nonlinearity, 4:961–979.
  * Sauer et al. (1997) Sauer, T., C. Grebogi, and J. A. Yorke (1997), How Long Do Numerical Chaotic Solutions Remain Valid?, Phys. Rev. Lett., 79:59–62, doi:10.1103/PhysRevLett.79.59.
  * Sauer (2002) Sauer, T. D. (2002), Shadowing breakdown and large errors in dynamical simulations of physical systems, Phys. Rev. E, 65:036,220, doi:10.1103/PhysRevE.65.036220.
  * Smith (2000) Smith, L. A. (2000), Disentangling uncertainty and error: On the predictability of nonlinear systems, in Nonlinear Dyanmics and Statistics, edited by A. Mees, pp. 31–64, Birkhäuser Boston.
  * Smith (2002) Smith, L. A. (2002), What might we learn from climate forecasts?, P. Natl Acad. Sci. USA, 99:2487–2492, doi:10.1073/pnas.012580599.
  * Smith and Spiegel (1987) Smith, L. A., and E. A. Spiegel (1987), Strange Accumulators, Ann. NY Acad. Sci., 497:61–65, doi:10.1111/j.1749-6632.1987.tb48711.x.
  * Smith et al. (1999) Smith, L. A., C. Ziehmann, and K. Fraedrich (1999), Uncertainty dynamics and predictability in chaotic systems, Q. J. Roy. Meteor. Soc., 125:2855–2886, doi:10.1002/qj.49712556005.
  * Smith et al. (2010) Smith, L. A., M. C. Cuéllar, H. Du, and K. Judd (2010), Exploiting Dynamical Coherence: A geometric approach to parameter estimation in nonlinear models, Phys. Lett. A, 374:2618–2623, doi:10.1016/j.physleta.2010.04.032.
  * Stemler and Judd (2009) Stemler, T., and K. Judd (2009), A guide to using shadowing filters for forecasting and state estimation, Physica D, 238:1260–1273, doi:10.1016/j.physd.2009.04.008.
  * Strogatz (1994) Strogatz, S. H. (1994), Nonlinear Dynamics and Chaos, Studies in Nonlinearity, 1st ed., Perseus Books Publishing LLC.
  * Stuart and Ord (1994) Stuart, A., and J. K. Ord (1994), Kendall’s advanced theory of statistics. Volume 1: Distribution Theory, 6th ed., Edward Arnold.
  * Toth and Kalnay (1993) Toth, Z., and E. Kalnay (1993), Ensemble Forecasting at NMC: The Generation of Perturbations, B. Am. Meteorol. Soc., 74:2317–2330, doi:10.1175/1520-0477(1993)074<2317:EFANTG>2.0.CO;2.
  * Winston and Krueger (1961) Winston, J. S., and A. F. Krueger (1961), Some Aspects of the Cycle of Available Potential Energy, Mon. Weather Rev., 89:307–318, doi:10.1175/1520-0493-89.9.307.
  * Wolf et al. (1985) Wolf, A., J. B. Swift, H. L. Swinney, and J. A. Vastano (1985), Determining Lyapunov exponents from a time series, Physica D, 16:285–317, doi:10.1016/0167-2789(85)90011-9.
  * Young and Read (2008) Young, R. M. B., and P. L. Read (2008), Flow transitions resembling bifurcations of the logistic map in simulations of the baroclinic rotating annulus, Physica D, 237:2251–2262, doi:10.1016/j.physd.2008.02.014.


canonical mapping functions, briefly as follows:

  * •

The two problems for $c_{1}(z)$ have the same extremal function $f(z)$. Such a
function maps $\Omega$ onto a circular slit disk centered at the origin and
having radius $\gamma_{\rm hydro}(z)$. More precisely, $f$ is given by

$f(z)=\exp\big{(}\gamma_{\rm hydro}(a)-{G}_{\rm
hydro}(z,a)-\mathrm{i}{G}^{*}_{\rm hydro}(z,a)\big{)}.$

  * •

The extremal function for the problem for $c_{D}(z)$ is a normalized version
of the integral of the hydrodynamic Bergman kernel:

$f(z)=\int_{a}^{z}\frac{{K_{\rm hydro}}(\zeta,a)}{{K_{\rm
hydro}}(a,a)}\,d\zeta.$

After mixing $K_{\rm hydro}$ (in suitable linear combinations) with the
corresponding adjoint kernel $L_{\rm hydro}$ one obtains canonical maps onto
horizontal and vertical slit domain.

  * •

The extremal function for the first problem for $c_{B}(z)$ is a suitably
normalized Ahlfors map:

$f(z)=\frac{K_{\text{Szeg\"{o}}}(z,a)}{K_{\text{Szeg\"{o}}}(a,a)L_{\text{Szeg\"{o}}}(z,a)}.$

This maps $\Omega$ onto the unit disk covered $\texttt{g}+1$ times. Here
$L_{\text{Szeg\"{o}}}(z,a)$ denotes the continuation of
$K_{\text{Szeg\"{o}}}(z,a)$ to the backside of the Schottky double as a half-
order differential. It has a simple pole at $z=a$ and is characterized by the
boundary relation

$K_{\text{Szeg\"{o}}}(z,a)\sqrt{dz}=\overline{L_{\text{Szeg\"{o}}}(z,a)}\sqrt{d\bar{z}}\quad{\rm
along}\,\,\partial\Omega.$

This is to be interpreted so that $\sqrt{dz}/\sqrt{d\bar{z}}=dz/ds$, the unit
tangent vector on $\partial\Omega$.

  * •

The extremal function for the second problem for $c_{B}(z)$ is the
(normalized) Szegö kernel:

$f(z)=\frac{K_{\text{Szeg\"{o}}}(z,a)}{K_{\text{Szeg\"{o}}}(a,a)}$

  * •

For $c_{\beta}$, the extremal function (multiple-valued) is the exponential of
the analytic completion of the Green function. More precisely:

$f(z)=\exp\big{(}\gamma_{\rm electro}(a)-G_{\rm electro}(z,a)-\mathrm{i}G_{\rm
electro}^{*}(z,a)\big{)}.$

  * •

For $M(z)$, the extremal function is the normalized Bergman kernel:

$f(z)=\frac{K_{\rm electro}(z,a)}{K_{\rm electro}(a,a)}.$

References for the inequalities (5.23), (5.24) are as follows.

  1. 1.

The first inequality is proved in [67] (p.172-179). It is actually obvious in
view of the extremal problems because a univalent function is by definition
analytic.

  2. 2.

The second inequality (related to Ahlfors-Beurling [1]) also appears in [67],
and in [65] (see also [66]) it is shown to be strict in the multiply connected
case (in [67] this is stated as an open question).

  3. 3.

The third inequality, in the form $\leq$, is obvious from the extremal
problems. Also the strict version follows easily, and is proved explicitly in
[67], p.178-179.

  4. 4.

The fourth inequality was conjectured in [72] and proved in [16]. Suita [72]
also has a simple proof of the inequality $4\pi
K_{\text{Szeg\"{o}}}(a,a)^{2}\leq K_{\rm electro}(a,a)$. Hejhal [48] shows
strictness inequality in the multiply connected case.

## References

  * [1] L. Ahlfors and A. Beurling, Conformal invariants and function-theoretic null-sets, Acta Math., 83 (1950), pp. 101–129.
  * [2] L. V. Ahlfors, Finitely generated Kleinian groups, Amer. J. Math., 86 (1964), pp. 413–429.
  * [3] , Conformal invariants: topics in geometric function theory, McGraw-Hill Series in Higher Mathematics, McGraw-Hill Book Co., New York-Düsseldorf-Johannesburg, 1973.
  * [4] L. V. Ahlfors and L. Sario, Riemann surfaces, Princeton Mathematical Series, No. 26, Princeton University Press, Princeton, N.J., 1960\.
  * [5] N. L. Alling and N. Greenleaf, Foundations of the theory of Klein surfaces, Lecture Notes in Mathematics, Vol. 219, Springer-Verlag, Berlin-New York, 1971.
  * [6] D. H. Armitage and S. J. Gardiner, Classical potential theory, Springer Monographs in Mathematics, Springer-Verlag London, Ltd., London, 2001\.
  * [7] V. I. Arnold, Mathematical methods of classical mechanics, Springer-Verlag, New York, 1978.  Translated from the Russian by K. Vogtmann and A. Weinstein, Graduate Texts in Mathematics, 60.
  * [8] V. I. Arnold and B. A. Khesin, Topological methods in hydrodynamics, vol. 125 of Applied Mathematical Sciences, Springer-Verlag, New York, 1998.
  * [9] M. F. Atiyah, Riemann surfaces and spin structures, Ann. Sci. École Norm. Sup. (4), 4 (1971), pp. 47–62.
  * [10] N. Balabanova and J. Montaldi, A hamiltonian approach for point vortices on non-orientable surfaces i: the möbius band, arXiv:2022.06160 [math.DS], 3 (2022), pp. 1–31.
  * [11] C. Bandle and M. Flucher, Harmonic radius and concentration of energy; hyperbolic radius and Liouville’s equations $\Delta U=e^{U}$ and $\Delta U=U^{(n+2)/(n-2)}$, SIAM Rev., 38 (1996), pp. 191–238.
  * [12] S. R. Bell, The Cauchy transform, potential theory and conformal mapping, Chapman & Hall/CRC, Boca Raton, FL, second ed., 2016.
  * [13] S. R. Bell and B. Gustafsson, Ruminations on Hejhal’s theorem about the Bergman and Szegö kernels, Anal. Math. Phys., 12 (2022), pp. Paper No. 24, 15.
  * [14] S. Bergman, The kernel function and conformal mapping, Mathematical Surveys, No. V, American Mathematical Society, Providence, R.I., revised ed., 1970\.
  * [15] S. Bergman and M. Schiffer, Kernel functions and conformal mapping, Compositio Math., 8 (1951), pp. 205–249.
  * [16] B. Berndtsson and L. Lempert, A proof of the Ohsawa-Takegoshi theorem with sharp estimates, J. Math. Soc. Japan, 68 (2016), pp. 1461–1472.
  * [17] S. Boatto and J. Koiller, Vortices on closed surfaces, in Geometry, mechanics, and dynamics, vol. 73 of Fields Inst. Commun., Springer, New York, 2015, pp. 185–237.
  * [18] L. A. Caffarelli and A. Friedman, Convexity of solutions of semilinear elliptic equations, Duke Math. J., 52 (1985), pp. 431–456.
  * [19] H. Cohn, Conformal mapping on Riemann surfaces, Dover Publications, Inc., New York, 1980.  Reprint of the 1967 edition, Dover Books on Advanced Mathematics.
  * [20] P. J. Davis, The Schwarz function and its applications, The Mathematical Association of America, Buffalo, N. Y., 1974.  The Carus Mathematical Monographs, No. 17.
  * [21] G. de Rham, Differentiable manifolds, vol. 266 of Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences], Springer-Verlag, Berlin, 1984.  Forms, currents, harmonic forms, Translated from the French by F. R. Smith, With an introduction by S. S. Chern.
  * [22] J. L. Doob, Classical potential theory and its probabilistic counterpart, vol. 262 of Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences], Springer-Verlag, New York, 1984\.
  * [23] P. L. Duren, Theory of $H^{p}$ spaces, Pure and Applied Mathematics, Vol. 38, Academic Press, New York-London, 1970.
  * [24] P. L. Duren and M. M. Schiffer, Robin functions and energy functionals of multiply connected domains, Pacific J. Math., 148 (1991), pp. 251–273.
  * [25] H. M. Farkas and I. Kra, Riemann surfaces, vol. 71 of Graduate Texts in Mathematics, Springer-Verlag, New York, second ed., 1992.
  * [26] H. Federer, Geometric measure theory, Die Grundlehren der mathematischen Wissenschaften, Band 153, Springer-Verlag New York Inc., New York, 1969.
  * [27] M. Fekete, Über die Verteilung der Wurzeln bei gewissen algebraischen Gleichungen mit ganzzahligen Koeffizienten, Math. Z., 17 (1923), pp. 228–249.
  * [28] M. Flucher, Variational problems with concentration, vol. 36 of Progress in Nonlinear Differential Equations and their Applications, Birkhäuser Verlag, Basel, 1999.
  * [29] M. Flucher and B. Gustafsson, Vortex motion in two-dimensional hydrodynamics, Royal Instituste of Technology Research Bulletins, TRITA-MAT-1997-MA-02 (1979), pp. 1–24.
  * [30] T. Frankel, The Geometry of Physics, Cambridge University Press, Cambridge, third ed., 2012.  An introduction.
  * [31] O. Frostman, Potentiel d’équilibre et capacité des ensembles avec quelques applications à la théorie des fonctions, Lunds Universitet, 1935.  Doctoral thesis.
  * [32] P. R. Garabedian, Schwarz’s lemma and the Szegö kernel function, Trans. Amer. Math. Soc., 67 (1949), pp. 1–35.
  * [33] P. Griffiths and J. Harris, Principles of algebraic geometry, Wiley-Interscience [John Wiley & Sons], New York, 1978.  Pure and Applied Mathematics.
  * [34] C. Grotta-Ragazzo, Errata and addenda to: “Hydrodynamic vortex on surfaces” and “The motion of a vortex on a closed surface of constant negative curvature” [3724607], J. Nonlinear Sci., 32 (2022), pp. Paper No. 63, 10.
  * [35] , Vortex on Surfaces and Brownian Motion in Higher Dimensions: Special Metrics, J. Nonlinear Sci., 34 (2024), p. Paper no. 31.
  * [36] C. Grotta-Ragazzo, B. Gustafsson, and J. Koiller, On the interplay between vortices and harmonic flows: Hodge decomposition of Euler’s equations in 2d, Regul. Chaotic Dyn., 29 (2024), pp. 241–303.
  * [37] C. Grotta Ragazzo and H. H. d. B. Viglioni, Hydrodynamic vortex on surfaces, J. Nonlinear Sci., 27 (2017), pp. 1609–1640.
  * [38] B. Gustafsson, On the motion of a vortex in two-dimensional flow of an ideal fluid in simply and multiply connected domains, Royal Instituste of Technology Research Bulletins, TRITA-MAT-1979-7 (1979), pp. 1–109.
  * [39] B. Gustafsson, Application of half-order differentials on Riemann surfaces to quadrature identities for arc-length, J. Analyse Math., 49 (1987), pp. 54–89.
  * [40] , On the convexity of a solution of Liouville’s equation, Duke Math. J., 60 (1990), pp. 303–311.
  * [41] , Lectures on balayage, in Clifford algebras and potential theory, vol. 7 of Univ. Joensuu Dept. Math. Rep. Ser., Univ. Joensuu, Joensuu, 2004, pp. 17–63.
  * [42] B. Gustafsson, Vortex motion and geometric function theory: the role of connections, Philos. Trans. Roy. Soc. A, 377 (2019), pp. 20180341, 27.
  * [43] , Vortex pairs and dipoles on closed surfaces, J. Nonlinear Sci., 32 (2022), pp. Paper No. 62, 38\.
  * [44] B. Gustafsson and J. Roos, Partial balayage on Riemannian manifolds, J. Math. Pures Appl. (9), 118 (2018), pp. 82–127.
  * [45] B. Gustafsson and A. Sebbar, Critical points of Green’s function and geometric function theory, Indiana Univ. Math. J., 61 (2012), pp. 939–1017.
  * [46] N. S. Hawley and M. Schiffer, Half-order differentials on Riemann surfaces, Acta Math., 115 (1966), pp. 199–236.
  * [47] H. Hedenmalm, B. Korenblum, and K. Zhu, Theory of Bergman spaces, vol. 199 of Graduate Texts in Mathematics, Springer-Verlag, New York, 2000.
  * [48] D. A. Hejhal, Theta functions, kernel functions, and Abelian integrals, American Mathematical Society, Providence, R.I., 1972.  Memoirs of the American Mathematical Society, No. 129.
  * [49] L. L. Helms, Potential theory, Universitext, Springer, London, second ed., 2014.
  * [50] J. Jost, Bosonic strings: a mathematical treatment, vol. 21 of AMS/IP Studies in Advanced Mathematics, American Mathematical Society, Providence, RI; International Press, Somerville, MA, 2001.
  * [51] B. Kawohl, Rearrangements and convexity of level sets in PDE, vol. 1150 of Lecture Notes in Mathematics, Springer-Verlag, Berlin, 1985.
  * [52] D. Kinderlehrer and G. Stampacchia, An introduction to variational inequalities and their applications, vol. 88 of Pure and Applied Mathematics, Academic Press, Inc. [Harcourt Brace Jovanovich, Publishers], New York-London, 1980.
  * [53] P. Koebe, Abhandlungen zur Theorie der konformen Abbildung, Acta Math., 41 (1916), pp. 305–344.  IV. Abbildung mehrfach zusammenhängender schlichter Bereiche auf Schlitzbereiche.
  * [54] , Abhandlungen zur Theorie der konformen Abbildung, Math. Z., 2 (1918), pp. 198–236.
  * [55] J. Koiller and S. Boatto, Vortex paires on surfaces, AIP Conference Proceedings 1130; https://doi.org/10.1063/1.3146241, 77 (2009).
  * [56] J. Koiller, B. Gustafsson, and C. Grotta-Ragazzo, On the interplay between vorticies and harmonic flows: Hodge decomposition of euler’s equations in 2d, arXiv:2309.12582, 77 (2023).
  * [57] C. C. Lin, On the Motion of Vortices in Two Dimensions, University of Toronto Studies, Applied Mathematics Series, no. 5, University of Toronto Press, Toronto, Ont., 1943.
  * [58] D. Minda and D. J. Wright, Univalence criteria and the hyperbolic metric, Rocky Mountain J. Math., 12 (1982), pp. 471–479.
  * [59] J.-J. Moreau, Fonctions convexes duales et points proximaux dans un espace hilbertien, C. R. Acad. Sci. Paris, 255 (1962), pp. 2897–2899.
  * [60] Z. Nehari, Conformal mapping, McGraw-Hill Book Co., Inc., New York, Toronto, London, 1952.
  * [61] K. Okikiolu, A negative mass theorem for surfaces of positive genus, Comm. Math. Phys., 290 (2009), pp. 1025–1031.
  * [62] A. K. Raina, Fay’s trisecant identity and conformal field theory, Comm. Math. Phys., 122 (1989), pp. 625–641.
  * [63] T. Ransford, Potential theory in the complex plane, vol. 28 of London Mathematical Society Student Texts, Cambridge University Press, Cambridge, 1995.
  * [64] E. B. Saff and V. Totik, Logarithmic potentials with external fields, vol. 316 of Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences], Springer-Verlag, Berlin, 1997\.  Appendix B by Thomas Bloom.
  * [65] M. Sakai, On constants in extremal problems of analytic functions, Kōdai Math. Sem. Rep., 21 (1969), pp. 223–225.
  * [66] , Quadrature domains, vol. 934 of Lecture Notes in Mathematics, Springer-Verlag, Berlin, 1982.
  * [67] L. Sario and K. Oikawa, Capacity functions, Die Grundlehren der mathematischen Wissenschaften, Band 149, Springer-Verlag New York Inc., New York, 1969.
  * [68] M. Schiffer and D. C. Spencer, Functionals of finite Riemann surfaces, Princeton University Press, Princeton, N. J., 1954.
  * [69] F. Schottky, Ueber die conforme Abbildung mehrfach zusammenhängender ebener Flächen, J. Reine Angew. Math., 83 (1877), pp. 300–351.
  * [70] H. S. Shapiro, The Schwarz function and its generalization to higher dimensions, University of Arkansas Lecture Notes in the Mathematical Sciences, 9, John Wiley & Sons Inc., New York, 1992.  A Wiley-Interscience Publication.
  * [71] J. Steiner, A geometrical mass and its extremal properties for metrics on $S^{2}$, Duke Math. J., 129 (2005), pp. 63–86.
  * [72] N. Suita, Capacities and kernels on Riemann surfaces, Arch. Rational Mech. Anal., 46 (1972), pp. 212–217.
  * [73] F. Trèves, Basic linear partial differential equations, Academic Press [A subsidiary of Harcourt Brace Jovanovich, Publishers], New York-London, 1975.  Pure and Applied Mathematics, Vol. 62.
  * [74] M. Tsuji, Potential theory in modern function theory, Chelsea Publishing Co., New York, 1975.  Reprinting of the 1959 original.
  * [75] H. Weyl, Die Idee der Riemannschen Fläche, Vierte Auflage. Unveränderter Nachdruck der dritten, Vollständig umgearbeiteten Auflage, B. G. Teubner Verlagsgesellschaft, Stuttgart, 1964.