# Prespecified-time observer-based distributed control of battery energy
storage systems

Wu Yang, Shu-Ming Liang, Yan-Wu Wang, and Zhi-Wei Liu
This work is supported by the National Natural Science Foundation of China
under Grant 61903147 and the Interdisciplinary Scientific Research Foundation
of GuangXi University under Grant 2022JCC023. Wu Yang is with the Key
Laboratory of Disaster Prevention and Structural Safety of Ministry of
Education and School of Electrical Engineering, Guangxi University, Nanning,
530004, P.R. China. (e-mail: biwuyang@163.com). Shu-Ming Liang is with the
School of Electrical Engineering, Guangxi University, Nanning, 530004, P.R.
China. (e-mail: a1017261884@163.com). Yan-Wu Wang and Zhi-Wei Liu are with the
School of Artificial Intelligence and Automation, Huazhong University of
Science and Technology, Wuhan, 430074, China. (e-mail<EMAIL_ADDRESS>zwliu@hust.edu.cn).

###### Abstract

This paper studies the state-of-charge (SoC) balancing and the total
charging/discharging power tracking issues for battery energy storage systems
(BESSs) with multiple distributed heterogeneous battery units. Different from
the traditional cooperative control strategies based on the asymptotical or
finite-time distributed observers, two distributed prespecified-time observers
are proposed to estimate average battery units state and average desired
power, respectively, which can be determined in advance and independent of
initial states or control parameters. Finally, two simulation examples are
given to verify the effectiveness and superiority of the proposed control
strategy.

###### Index Terms:

Battery energy storage system; State-of-charge; Power tracking; Distributed
control; Prespecified-time observer.

##  I Introduction

Energy storage system is an indispensable part of microgrid, which can
guarantee the power quality and reliability and reduce the energy lose[1].
Among all kinds of energy storage technologies such as supercapacitor,
superconducting magnetic, etc, BESS has an irreplaceable position in microgrid
because of its advantages of fast response speed, high energy density, high
efficiency and flexible configuration[2, 3]. During past decades, various
types of BESSs are rapidly integrated into microgrid[4, 5]. Despite their
continuous advances in electrochemical technology, the control of BESSs
remains a very challenging problem[6].

Generally, BESSs are composed of multiple battery units, which can communicate
with its nearby battery units and monitor and control the charging/discharging
power of itself. The primary goal of BESS management is to balance the SoC of
all battery units and to meet the total charging/discharging desired power of
microgrid. Note that, due to variations and deviations in manufacturing
processes and operating conditions, battery units can manifest different
characteristics even with the same specifications. Consequently, how to design
an appropriate control strategies for BESS has attracted the broad interest of
researchers [7, 8, 9, 10, 11].

The control methods for BESS are typically divided into centralized,
decentralized or distributed ones: in centralized control method, a central
controller is needed to monitor and coordinate the SoC and other critical
states of all battery units[12], which is costly and prone to single-point
failure; Compared with centralized control method, neither central controller
nor communication between each battery units is required in decentralized
control method, however, it has been reported in [13] that such control method
may lead to slow SoC balancing, low SoC balancing accuracy and poor bus
voltage quality; By contrast, in distributed control method, no extra
centralized controller is needed and each battery unit can communicate with
its neighbors, such method not only can save communication cost, but also has
the advantages of robustness and reconfigurability [14]. Motivated by this,
different distributed cooperative control schemes are developed for BESSs. For
examples, [8] proposed a distributed SoC balancing control method based on
event-triggered mechanism, which allows each battery unit transmits signals to
its neighboring ones only when the triggered condition is met. Unfortunately,
the average desired power and average battery unit state are needed when
implementing the proposed controller, which is unrealistic since both the
average desired power and average battery unit state are of the global
information rather than the local one. To overcome drawback in [8], recently,
[11] introduced two types of distributed observers, one can achieve asymptotic
estimation on the average desired power and average battery unit state, while
the other ensures finite-time estimation. It should be emphasized that the
proposed distributed asymptotic observers in [11] can only be carried out as
time tends to infinity. In fact, finite-time convergence is more desirable in
applications due to the advantage of faster convergence rate and robustness
against uncertainties. Although the distributed finite-time observers are
further constructed in [11], the observer convergence rate depends on the
initial values and certain parameters of BESSs, which may restrict its
applications since the settling time is not fixed for different initial
values. Note that in the framework of prespecified-time stability, the
settling time can be bounded by a fixed value, independent of the initial
conditions and can be pre-set according to the task requirements, which makes
prespecified-time stability more desirable. To the best of our knowledge, the
SoC balancing and the total charging/discharging power tracking issues for
heterogeneous BESSs has not been fully investigated.

Inspired by the above observations, this paper addresses the problem of the
SoC balancing and the total charging/discharging power tracking for
heterogeneous BESS. The main contribution of this paper is highlighted as
follows: Different from the existing control strategies [8, 11], the
prespecified-time distributed observer-based cooperative control method is
employed to achieve both the SoC balancing and the total charging/discharging
power tracking for heterogeneous BESS. Specifically, by introducing time-
dependent function, two prespecified-time distributed observers are proposed
to reconstruct the average desired power and the average battery unit state,
respectively. It is shown that the settling time of the proposed observer can
be pre-set by user. Moreover, simulation examples are given to verify that the
performance of the proposed control method is better than that of [8, 11].

The rest of the paper is organized as follows. Section II reviews graph theory
and describes BESSs. Section III introduces the charging/discharging power
controller based on the prespecified-time average power observer and the
prespecified-time average battery unit state observer. Section IV provides
simulation results. Section V concludes this paper.

##  II Preliminaries and system description

###  II-A Preliminaries

In this subsection, graph theory is employed to describe the communication of
BESSs, in which vertex $i\in\mathcal{V}$ refers to battery unit
$i\in\mathcal{V}$ and edge $(i,j)\in\mathcal{E}$ refers to the interconnection
between the battery units $i$ and $j$, where $\mathcal{V}$ and $\mathcal{E}$
denote the vertex set and edge set, respectively. The adjacency matrix
$A=[a_{ij}]$ associated with graph $\mathcal{G}$ is defined as
$a_{ij}=a_{ji}=0$ if there is no communication link between unit $i$ and unit
$j$, otherwise $a_{ij}=a_{ji}=1$. The neighboring set of agent $i$ is
represented by
$\mathcal{N}_{i}=\left\\{j\in\mathcal{V}:(i,j)\in\mathcal{E}\right\\}$. The
Laplacian matrix of graph $\mathcal{G}$ is defined as $L=[l_{ij}]$, where
$l_{ij}=-a_{ij}$ if $i\neq j$ and $l_{ii}=\sum_{k=1,k\neq i}^{N}a_{ik}$.
Obviously, the Laplacian $L$ is symmetric with eigenvalues
$0=\lambda_{1}(L)<\lambda_{2}(L)\leq\lambda_{3}\leq\cdots\leq\lambda_{N}(L)$.
A graph is connected if there exists a path between any two distinct nodes.

###  II-B System description

Consider a class of heterogeneous BESS with $N$ battery units, described as
follows based on the Coulomb counting method,

$\begin{array}[]{l}s_{i}(t)=s_{i}(0)-\frac{1}{C_{i}}\int_{0}^{t}i_{i}(\tau)d\tau,\\\
p_{i}(t)=V_{i}(t)i_{i}(t),i\in\mathcal{N}=\\{1,2,\cdots,N\\},\end{array}$ (1)

where $s_{i}(t)$, $s_{i}(0)$, $C_{i}$ and $i_{i}(t)$ are the SoC value, the
initial SoC value, the capacity and the output current of the $i$-th battery
unit, respectively. $p_{i}(t)$ is the output power of the $i$-th battery unit,
$V_{i}(t)$ is the output voltage. Note that $p_{i}<0$ indicates changing and
$p_{i}>0$ indicates dischanging. In general, the output voltage $V_{i}(t)$ of
each unit is assumed to be unchanged during its operation[15]. Differentiating
both sides of the first equation of (1) yields

$\dot{s}_{i}(t)=-\frac{1}{C_{i}}i_{i}(t),i\in\mathcal{N}.$ (2)

Then, from (1) and (2), one has

$\dot{s}_{i}(t)=-\frac{1}{C_{i}V_{i}}p_{i}(t),i\in\mathcal{N}.$ (3)

For each battery unit $i\in\mathcal{N}$, define the battery unit state

$x_{i}(t)=\left\\{\begin{array}[]{l}C_{i}V_{i}s_{i}(t),\text{ discharging
mode},\\\ C_{i}V_{i}(1-s_{i}(t)),\text{ charging mode}.\end{array}\right.$ (4)

Moreover, the following assumption is made.

######  Assumption 1 ([11])

There exist constants $a_{1},a_{2}>0$ such that

$a_{1}\leq x_{i}(t)\leq a_{2},t\geq 0,i\in\mathcal{N}.$

###### Remark 1

Assumption 1 ensures that SoC of battery units varies within an appropriate
range to avoid overcharging or overdischarging during the operation of BESSs.

Then it deduces from (3) that

$\dot{x}_{i}(t)=\left\\{\begin{array}[]{l}C_{i}V_{i}\dot{s}_{i}(t)=-p_{i}(t),\\\
-C_{i}V_{i}\dot{s}_{i}(t)=p_{i}(t).\end{array}\right.$ (5)

Denote the average battery unit state and the average desired power as,
respectively,

$x_{a}(t)=\frac{1}{N}\sum_{i=1}^{N}x_{i}(t),\\\ p_{a}(t)=\frac{1}{N}p^{*}(t),$

where $p^{*}(t)$ is total desired power.

In this paper, the following assumption on $p^{*}(t)$ is required.

######  Assumption 2 ([8, 11])

The total desired charging/discharging power $p^{*}(t)$ of BESSs satisfies

$\underline{P}\leq\mid p^{*}(t)\mid\leq\bar{P},t\geq 0,\\\
\mid\dot{p}^{*}(t)\mid\leq\epsilon,t\geq 0,$

where $\bar{P}$, $\underline{P}$, and $\epsilon$ are positive constants.

In this paper, each battery unit is able to communicate with its neighbors,
then the following assumptions on the communication topology of BESSs are
made.

######  Assumption 3 ([8, 11])

The communication topology of BESSs is an undirected connected graph.

######  Assumption 4 ([11])

At least one battery unit has access to the total desired power $p^{*}(t)$.

From Assumptions 3 and 4, the matrix $H=L+B>0$[16], where the diagonal matrix
$B={\rm diag}\\{b_{1},b_{2},\cdots,b_{N}\\}$, $b_{i}=1$ if the $i$-th battery
unit gets access to the total desired power $p^{*}(t)$ and $b_{i}=0$
otherwise.

The control objective of BESS is stated as follows.

###### Problem 1

Consider BESSs (1), assume that Assumptions1-4 hold. By using the local
information, design the SoC balancing and the charging/discharging power
tracking distributed controllers for the battery units so that

  1. 1.

all battery units achieve SoC balancing with any pre-specified accuracy
$\varepsilon_{soc}\geq 0$ , i.e.,

$\lim_{t\to\infty}\|s_{i}(t)-s_{j}(t)\|\leq\varepsilon_{soc},i,j\in\\{1,2,\cdots,N\\};$

  2. 2.

the total charging/discharging power of BESSs tracks the total desired power
with any pre-specified accuracy $\varepsilon_{power}\geq 0$, i.e.,

$\lim_{t\to\infty}\|\sum\limits_{i=1}^{N}p_{i}-p^{*}(t)\|\leq\varepsilon_{power},i\in\\{1,2,\cdots,N\\}.$

##  III Power Control Based On Prespecified-time observer

To solve Problem 1, the following control strategy is developed in [8],

$p_{i}(t)=\frac{x_{i}(t)}{x_{a}(t)}p_{a}(t).$ (6)

Note that both $p_{a}(t)$ and $x_{a}(t)$ are global information, which are not
available to every battery unit. Recently, [11] introduced two types of
distributed observers, one can achieve asymptotic estimation on $p_{a}(t)$ and
$x_{a}(t)$, while the other ensures finite-time estimation. Unfortunately, the
proposed distributed asymptotic observers in [11] can only be carried out as
time tends to infinity, while the observer convergence rate of the distributed
finite-time observers depends on the initial values of BESSs. In fact, the
prespecified-time convergence property is more desirable. Motivated by [17],
firstly, the following prespecified-time distributed observer is proposed to
estimate $p_{a}(t)$,

$\begin{array}[]{l}\dot{\hat{p}}_{a,i}=-\alpha{\rm
sign}(v_{i})-\psi\frac{r}{t_{b}-t_{0}}\omega(t)v_{i},\\\
v_{i}=\sum_{j\in{N_{i}}}(\hat{p}_{a,i}-\hat{p}_{a,j})+b_{i}(\hat{p}_{a,i}-p_{a}),\end{array}$
(7)

where $\hat{p}_{a,i}$ is the $i$-th battery unit’s estimation of $p_{a}(t)$,
$\alpha,\psi$ and $r$ are positive parameters, $v_{i}$ is the internal state,
the function $\omega(t)$ is defined as follows:

$\omega(t)=\begin{cases}\frac{t_{b}-t_{0}}{t_{b}-t}&t\in\left[t_{0},t_{b}\right),\\\
1&t\in\left[t_{b},\infty\right),\end{cases}$

where $t_{b}$ is a prespecified time.

###### Lemma 1

Consider BESSs (1) with the prespecified-time distributed observer (7). Assume
that Assumptions 2-4 hold. Then, the average desired power $p_{a}(t)$ can be
estimated accurately at prespecified time $t_{b}$, i.e.,
$\hat{p}_{a,i}(t)=p_{a}(t)$ over $t\in[t_{b},+\infty)$, if the observer
parameter $\alpha$ is selected such that $\alpha\geq\frac{\epsilon}{N}$.

###### Proof 1

Denote $\hat{p}_{a}=(\hat{p}_{a,1},\hat{p}_{a,2},\cdots,\hat{p}_{a,N})^{T}$,
$\tilde{p}=\hat{p}_{a}-1_{N}p_{a},$ and
$\tilde{p}=(\tilde{p}_{1},\tilde{p}_{2},\cdots,\tilde{p}_{N})^{T}$, then it
follows from (7) that

$\displaystyle v$ $\displaystyle=L\hat{p}_{a}+B\hat{p}_{a}-B1_{N}p_{a}$
$\displaystyle=H\hat{p}_{a}-B1_{N}p_{a}$ $\displaystyle=H\tilde{p},$

where $H=L+B$.

Consider the Lyapunov function candidate $V_{1}(t)=\frac{1}{2}v^{T}H^{-1}v$.
Then, The time derivative of $V_{1}(t)$ is evaluated as

$\displaystyle\dot{V}_{1}(t)$ $\displaystyle=v^{T}H^{-1}H\dot{\tilde{p}}$
$\displaystyle=v^{T}(\dot{\hat{P}}_{a}-1_{N}\frac{1}{N}\dot{p}^{*})$
$\displaystyle\leq-\alpha\sum_{i=1}^{N}|v_{i}|-\frac{r\psi}{t_{b}-t_{0}}\omega(t)\sum_{i=1}^{N}v_{i}^{2}+\frac{\epsilon}{N}\|v\|_{1}$
$\displaystyle\leq(-\alpha+\frac{\epsilon}{N})\|v\|_{1}-\frac{r\psi}{t_{b}-t_{0}}\omega(t)\|v\|_{2}^{2}$
$\displaystyle\leq-\frac{r\psi}{t_{b}-t_{0}}\omega(t)\|v\|_{2}^{2}$
$\displaystyle\leq-\frac{r\psi}{t_{b}-t_{0}}\omega(t)\frac{2}{\lambda_{\max}(H^{-1})}V_{1}(t),$
(8)

by using the fact that $\alpha\geq\frac{\epsilon}{N}$ and
$\frac{1}{2}v^{T}H^{-1}v\leq\frac{1}{2}\lambda_{\max}(H^{-1})\|v\|_{2}^{2}.$

Denote $\Omega(t)=\omega^{r}(t)$, it gives that

$\frac{r}{t_{b}-t_{0}}\omega(t)=\frac{\dot{\Omega}(t)}{\Omega(t)}.$ (9)

According to (1) and (9), it gets that

$V_{1}(t)\leq[\omega(t)]^{-\frac{2\psi
r}{\lambda_{\max}(H^{-1})}}[\omega(t_{0})]^{\frac{2\psi
r}{\lambda_{\max}(H^{-1})}}V_{1}(t_{0}).$ (10)

Moreover, one from the definition of $\omega(t)$ that

$[\omega(t_{0})]^{\frac{2\psi r}{\lambda_{\max}(H^{-1})}}=1,\lim_{t\rightarrow
t^{-}_{b}}[\omega(t)]^{-\frac{2\psi r}{\lambda_{\max}(H^{-1})}}=0,$ (11)

which implies that

$0\leq\lim_{t\rightarrow t^{-}_{b}}V_{1}(t)\leq 0.$ (12)

Therefore, $\lim_{t\rightarrow t^{-}_{b}}\|v\|=0$, so does $\lim_{t\rightarrow
t^{-}_{b}}\|\tilde{p}\|=0$, that is, the average desired power $p_{a}(t)$ is
estimated accurately at prespecified time $t_{b}$.

Furthermore, for $t\in[t_{b},+\infty)$, following the similar step in above,
one obtains that

$\dot{V}_{1}(t)\leq-\frac{r\psi}{t_{b}-t_{0}}\|v\|_{2}^{2}.$

Since $v$ is continuous, $V_{1}(t)$ is also continuous. Therefore,
$\lim_{t\rightarrow t_{b}}V_{1}(t)=V_{1}(t_{b})$. Then one can get

$0\leq V_{1}(t)\leq V_{1}(t_{b})=0,t\in[t_{b},+\infty),$

which means that $V_{1}(t)\equiv 0$ on $t\in[t_{b},+\infty)$, thus
$\tilde{p}\equiv 0$ and $\dot{\hat{p}}_{a}\equiv 0$ over
$t\in[t_{b},+\infty)$.

Lastly, the boundedness of the proposed observer over $[t_{0},\infty)$ is
checked.

Denote $\tilde{\psi}_{1}=\frac{\psi}{\lambda_{\max}(H^{-1})}$, for
$t\in[t_{0},t_{b})$, it can be obtained from (10) that

$\|\tilde{p}(t)\|\leq\sqrt{\frac{\lambda_{\max}(H)}{\lambda_{\min}(H)}}[\omega(t)]^{-\tilde{\psi}_{1}r}\|\tilde{p}(t_{0})\|$
(13)

and

$\|\omega(t)\tilde{p}(t)\|\leq\sqrt{\frac{\lambda_{\max}(H)}{\lambda_{\min}(H)}}[\omega(t)]^{-(\tilde{\psi}_{1}r-1)}\|\tilde{p}(t_{0})\|.$
(14)

When $\tilde{\psi}_{1}r>2$, one has that
$0<[\omega(t)]^{-\tilde{\psi}_{1}r}<1$ and
$0<[\omega(t)]^{-(\tilde{\psi}_{1}r-1)}<1$. Then

$\max\\{\|\omega(t)\tilde{p}(t)\|,\|\tilde{p}(t)\|\\}\leq\
\sqrt{\frac{\lambda_{\max}(H)}{\lambda_{\min}(H)}}\|\tilde{p}(t_{0})\|.$ (15)

Moreover, it follows from (7) that

$\displaystyle\|\dot{\hat{p}}_{a}\|$ $\displaystyle\leq\
\alpha\|H\|\|\tilde{p}\|+\frac{r\psi}{t_{b}-t_{0}}\|H\|\|\omega(t)\tilde{p}\|$
$\displaystyle\leq(\alpha+\frac{r\psi}{t_{b}-t_{0}})\|H\|\sqrt{\frac{\lambda_{\max}(H)}{\lambda_{\min}(H)}}\|\tilde{p}_{0}\|,$
(16)

which confirms that $\hat{p}_{a}$ is bounded on $[t_{0},t_{b})$. Similarly, it
can be seen from the definition that $\hat{p}_{a}(t)$ is bounded on
$[t_{b},\infty)$. $\blacksquare$

For each battery unit $i\in\mathcal{N}$, to estimate $x_{a}(t)$, the following
prespecified-time battery unit average state observer is further constructed:

$\displaystyle\dot{q}_{i}=$ $\displaystyle-\beta{\rm
sign}(\hat{\xi}_{i})-\psi\frac{r}{t_{a}}\omega(t)\hat{\xi}_{i}$
$\displaystyle\hat{\xi}_{i}=$
$\displaystyle\sum_{j=1}^{n}a_{ij}(\hat{x}_{a,i}-\hat{x}_{a,j})$
$\displaystyle\hat{x}_{a,i}=$
$\displaystyle\sum_{j=1}^{n}a_{ij}(q_{i}-q_{j})+x_{i}$ (17)

where $\hat{x}_{a,i}$ is the estimate state of the $i$-th battery unit for
$x_{a}(t)$, $\beta$ is a positive parameter, $\hat{\xi}_{i}$ and $q_{i}$ are
the internal states.

###### Lemma 2

Consider BESSs (1) with the prespecified-time distributed observer (III).
Assume that Assumptions 1, 3 and 4 hold. Then, the average battery unit state
$x_{a}(t)$ can be estimated accurately at prespecified time $t_{b}$, i.e.,
$\hat{x}_{a,i}(t)=x_{a}(t)$ over $t\in[t_{b},+\infty)$, if the observer
parameter $\beta$ is selected such that
$\beta\geq\frac{\bar{P}}{\sqrt{N}\lambda_{2}(L)}$.

###### Proof 2

Denote $\hat{x}_{a}=(\hat{x}_{a,1},\hat{x}_{a,2},\cdots,\hat{x}_{a,N})^{T}$,
$\tilde{x}=\hat{x}_{a}-x_{a}1_{N}$ and
$\tilde{x}=(\tilde{x}_{1},\tilde{x}_{2},\cdots,\tilde{x}_{N})^{T}$, then one
can obtain that

$\displaystyle\dot{\tilde{x}}$
$\displaystyle=\dot{\hat{x}}_{a}-\dot{x}_{a}1_{N}$
$\displaystyle=L\dot{q}+\dot{x}-\dot{x}_{a}1_{N}$
$\displaystyle=L\dot{q}+I_{n}\dot{x}-\frac{1}{n}1_{N}1_{N}^{T}\dot{x}$
$\displaystyle=L\dot{q}+LL^{\dagger}\dot{x},$ (18)

where $L^{\dagger}$ is the generalized inverse matrix of $L$.

Consider the Lyapunov function candidate
$V_{2}(t)=\frac{1}{2}\tilde{x}^{T}\tilde{x}$. The derivative of $V_{2}(t)$
along the trajectory of (2) is evaluated as

$\displaystyle\dot{V}_{2}(t)$
$\displaystyle=\tilde{x}^{T}\dot{\tilde{x}}=\tilde{x}^{T}L\dot{q}+\tilde{x}^{T}LL^{\dagger}\dot{x}$
$\displaystyle=\mu^{T}(t)\dot{q}+\mu^{T}(t)L^{\dagger}\dot{x},$ (19)

where $\mu(t)=L\tilde{x}$.

In fact,

$\displaystyle\mu^{T}(t)\dot{q}$
$\displaystyle=-\sum_{i=1}^{n}\quad\beta|\mu_{i}(t)|-\mu^{T}(t)\psi\frac{r}{t_{b}-t_{0}}\omega(t)L\hat{x}_{a}$
$\displaystyle\leq-\beta\|\mu(t)\|_{1}-\psi\frac{r}{t_{b}-t_{0}}\omega(t)\tilde{x}^{T}L^{T}L(\tilde{x}+x_{a}1_{N})$
$\displaystyle\leq-\beta\|\mu(t)\|_{1}-\psi\frac{r}{t_{b}-t_{0}}\omega(t)\tilde{x}^{T}L^{T}L\tilde{x}.$
(20)

and

$\displaystyle\mu^{T}(t)L^{\dagger}\dot{x}$
$\displaystyle\leq\sigma_{\max}(L^{\dagger})\|\mu(t)\|_{2}\|\dot{x}(t)\|_{2}$
$\displaystyle\leq\sigma_{\max}(L^{\dagger})\sqrt{N}\|\mu(t)\|_{1}\|\dot{x}(t)\|_{\infty}$
$\displaystyle\leq\frac{\sqrt{N}}{\lambda_{2}(L)}\|\mu(t)\|_{1}\|\dot{x}(t)\|_{\infty}$
$\displaystyle\leq\frac{\sqrt{N}\bar{P}}{N\lambda_{2}(L)}\|\mu(t)\|_{1},$ (21)

where $\sigma_{\max}(\cdot)$ denotes the maximum singular value. Then, one can
further get

$\displaystyle\dot{V}_{2}(t)$
$\displaystyle\leq(\frac{\sqrt{N}\bar{P}}{N\lambda_{2}(L)}-\beta)\|\mu(t)\|_{1}-\psi\frac{r}{t_{a}}\omega(t)\tilde{x}^{T}L^{T}L\tilde{x}$
$\displaystyle\leq-\psi\frac{r}{t_{a}}\omega(t)\tilde{x}^{T}L^{T}L\tilde{x},$
(22)

by using $\beta\geq\frac{\bar{P}}{\sqrt{N}\lambda_{2}(L)}$.

Since the matrix $L$ is a symmetrical, then there exists a matrix $T$ such
that $L=T^{H}\Lambda T$. Let $Y=T^{H}\tilde{x}$, one can deduce that

$\displaystyle\dot{V}_{2}(t)$
$\displaystyle\leq-\psi\frac{r}{t_{b}-t_{0}}\omega(t)Y^{T}T^{H}L^{T}LTY$
$\displaystyle\leq-\psi\frac{r}{t_{b}-t_{0}}\omega(t)\lambda_{2}(L^{T}L)Y^{T}Y$
$\displaystyle\leq-2\psi\frac{r}{t_{b}-t_{0}}\omega(t)\lambda_{2}(L^{T}L)V_{2}(t)$
$\displaystyle=-2\psi\frac{\dot{\Omega}(t)}{\Omega(t)}\lambda_{2}(L^{T}L)V_{2}(t).$
(23)

Which means that

$V_{2}(t)\leq[\omega(t)]^{-2r\psi\lambda_{2}(L^{T}L)}[\omega(t_{0})]^{2r\psi\lambda_{2}(L^{T}L)}V_{2}(t_{0}),$
(24)

and

$0\leq\lim_{t\rightarrow t^{-}_{b}}V_{2}(t)\leq 0.$ (25)

Therefore, $\lim_{t\rightarrow t^{-}_{b}}\|\tilde{x}\|=0$, that is, the
average battery unit state $x_{a}(t)$ is estimated accurately at prespecified
time $t_{b}$. Following the similar step in Lemma 1, for
$t\in[t_{b},+\infty)$, one obtains that $V_{2}(t)\equiv 0$ on
$t\in[t_{b},+\infty)$, thus $\tilde{x}\equiv 0$ over $t\in[t_{b},+\infty)$.

Lastly, the boundedness of $\hat{x}_{a,i}$ over $[t_{0},\infty)$ is verified.

Denote $\tilde{\psi}_{2}=\psi\lambda_{2}(L^{T}L)$, for $t\in[t_{0},t_{b})$, it
can be obtained from (24) that

$\|\tilde{x}(t)\|\leq[\omega(t)]^{-\tilde{\psi}_{2}r}\|\tilde{x}(t_{0})\|$
(26)

and

$\|\omega(t)\tilde{x}(t)\|\leq[\omega(t)]^{-(\tilde{\psi}_{2}r-1)}\|\tilde{x}(t_{0})\|.$
(27)

Moreover, it can be derived that

$\max\\{\|\omega(t)\tilde{x}(t),\|\tilde{x}(t)\|\\}\|\leq\
\|\tilde{x}(t_{0})\|$ (28)

when $\tilde{\psi}_{2}r>2$.

Therefore,
$\|\dot{q}\|\leq\beta+n\frac{r\psi}{t_{a}}\|L\|\|\tilde{x}(t_{0})\|,$ which
means that $q$ is bounded on $[t_{0},t_{b})$, and $\|\hat{x}_{a}\|\leq\
n\|L\|\|q\|+n\|x\|,$ which confirms that $\hat{x}_{a}$ is bounded on
$[t_{0},t_{b})$. Similarly, it can be seen from the definition that
$\hat{x}_{a}(t)$ is bounded on $[t_{b},\infty)$. $\blacksquare$

Based on Lemmas 1 and 2, the following result on the performance of the power
control algorithms (6) with the prespecified-time average power observer (7)
and the prespecified-time battery unit state observer (III) is given.

###### Theorem 1

Suppose that Assumptions 1-4 hold for BESSs (1), the observer parameters
$\alpha$ and $\beta$ are selected as in Lemmas 1 and 2, respectively, then,
under the power control algorithms (6) with the prespecified-time distributed
observers (7) and (III), it can be concluded that Problem 1 is solved.

###### Proof 3

Based on Lemmas 1 and 2, $\hat{p}_{a,i}(t)=p_{a}(t)$,
$\hat{x}_{a,i}(t)=x_{a}(t)$, when $t\in[t_{b},+\infty)$. Then, similar to [8],
both the SoC balancing and the desired power tracking are achieved by the
controller (6). $\blacksquare$

###### Remark 2

In contrast to the asymptotic observers-based cooperative controller or the
finite-time observers-based cooperative controller[11], the proposed
cooperative controller relies on the prespecified-time observer, where the
time-dependent function $\omega(t)$ is leveraged.

###### Remark 3

The function $\omega(t)$ plays a key role in ensuring the prespecified-time
stability of the corresponding estimation error systems. Specifically, from
the definition of $\omega(t)$, one has that $\lim_{t\rightarrow
t^{-}_{b}}[\omega(t)]^{-\tau}=0$ and $[\omega(t_{0})]^{\tau}=1$ for any
$\tau>0$, which leads to $\lim_{t\rightarrow t^{-}_{b}}V_{i}(t)=0,i=1,2$.

###### Remark 4

Compared with the asymptotic observers and the finite-time observers presented
in [11], the converge property of the former can only be ensured when
$t\rightarrow\infty$, while the converge rate of the latter depends on the
initial state of BESSs, the prespecified time observers are designed herein,
and its converge rate can be determined in advance according to task
requirement. Moreover, in the next section, some simulation examples are given
to show the advantage of our results.

##  IV Simulation Examples

In this section, the effectiveness and superiority of the proposed control
strategy with some existing works, such as [11], are demonstrated by numerical
examples.

###### Example 1

Consider the following BESSs with $6$ battery units, the parameters of the
battery units are $V=20$V, $C=220$Ah, and the initial SoC of the battery units
are $(0.95,0.86,0.83,0.93,0.97,0.88)$. The communication topology of BESS is
shown in Fig. 1. Besides, assume that only battery unit $1$ has access to the
desired charging/discharging power, i.e., $b_{1}=1$, and
$b_{i}=0,i=2,3,\cdots,6$. Moreover, assume that the desired power is given as
follows:

  1. 1.

Case 1: $p^{*}(t)=\pm(4200\sin(t)+4200$)W.

  2. 2.

Case 2:

$p^{*}(t)=\pm\begin{cases}1000&\text{ $t\in[0,0.25)$},\\\
1500\sin(5t)+2000&\text{ $t\in[0.25,0.5)$},\\\ 5000t&\text{
$t\in[0.5,0.75)$},\\\ -4000t+6000&\text{ $t\in[0.75,1)$}.\end{cases}$

Figure 1: The communication network topology

For the prespecified-time distributed observers (7) and (III), the
prespecified time is chosen as $t_{b}=4$ and $t_{0}=0$, the observer
parameters are selected as $\psi=4$, $r=50$, $\alpha=1000$, $\beta=3430$, the
initial states of the observers (7) and (III) are given as $v_{i}(0)=0$,
$\hat{p}_{i}(0)=p_{i}(0)=0$, $q_{i}(0)=0$, $\hat{x}_{i}(0)=x_{i}(0)$,
$i=1,2,\cdots,6$. Then, the simulation results under the controller (6) with
the observers (7) and (III) are shown in Figs.2-4. Where the solid line
represents the simulation results under our control method, while the dotted
line represents that in [11]. It can be observed that the proposed cooperative
control method is better in comparison with the scheme in [11].

(a) SoC of all battery units. (b) Output power of all battery units.

Figure 2:  BESSs under the controller (6) with the observers (7) and (III) [In
discharge mode, Case 1]

(a) SoC of all battery units. (b) Output power of all battery units.

Figure 3:  BESSs under the controller (6) with the observers (7) and (III) [In
charge mode, Case 1]

(a) Output power of all battery units (In discharge mode). (b) Output power of
all battery units (In charge mode).

Figure 4:  BESSs under controllers (6) with the observers (7) and (III) [Case
2]

##  V Conclusions

This paper has investigated the problem of SoC balancing and power tracking
control for BESSs with multiple heterogenous battery units. Two prespecified-
time distributed observers have been constructed to estimate the average
desired power and the average battery unit state, respectively. Then, by
selecting observer’s parameters appropriately, the effectiveness of the
proposed control strategy based on the prespecified-time distributed observers
has been further analyzed and has been verified by some simulation examples.
Future works include extending this results to BESSs with time delays or more
communication constraints.

## References

  * [1] H. Chen, T. N. Cong, W. Yang, C. Tan, Y. Li, and Y. Ding, “Progress in electrical energy storage system: A critical review,” Progress in Natural Science, vol. 19, no. 3, pp. 291–312, 2009.
  * [2] M. T. Lawder, B. Suthar, P. W. Northrop, S. De, C. M. Hoff, O. Leitermann, M. L. Crow, S. Santhanagopalan, and V. R. Subramanian, “Battery energy storage system (bess) and battery management system (bms) for grid-scale applications,” Proceedings of the IEEE, vol. 102, no. 6, pp. 1014–1030, 2014.
  * [3] C. A. Hill, M. C. Such, D. Chen, J. Gonzalez, and W. M. Grady, “Battery energy storage for enabling integration of distributed solar power generation,” IEEE Transactions on Smart Grid, vol. 3, no. 2, pp. 850–857, 2012.
  * [4] T. Feehally, A. Forsyth, R. Todd, M. Foster, D. Gladwin, D. Stone, and D. Strickland, “Battery energy storage systems for the electricity grid: Uk research facilities,” 2016.
  * [5] B. M. Gundogdu, S. Nejad, D. T. Gladwin, M. P. Foster, and D. A. Stone, “A battery energy management strategy for uk enhanced frequency response and triad avoidance,” IEEE Transactions on Industrial Electronics, vol. 65, no. 12, pp. 9509–9517, 2018.
  * [6] H. Rahimi-Eichi, U. Ojha, F. Baronti, and M.-Y. Chow, “Battery management system: An overview of its application in the smart grid and electric vehicles,” IEEE Industrial Electronics Magazine, vol. 7, no. 2, pp. 4–16, 2013.
  * [7] H. Cai and G. Hu, “Distributed control scheme for package-level state-of-charge balancing of grid-connected battery energy storage system,” IEEE Transactions on industrial informatics, vol. 12, no. 5, pp. 1919–1929, 2016.
  * [8] L. Xing, Y. Mishra, Y.-C. Tian, G. Ledwich, C. Zhou, W. Du, and F. Qian, “Distributed state-of-charge balance control with event-triggered signal transmissions for multiple energy storage systems in smart grid,” IEEE Transactions on Systems, Man, and Cybernetics: Systems, vol. 49, no. 8, pp. 1601–1611, 2019.
  * [9] C. Deng, Y. Wang, C. Wen, Y. Xu, and P. Lin, “Distributed resilient control for energy storage systems in cyber–physical microgrids,” IEEE Transactions on Industrial Informatics, vol. 17, no. 2, pp. 1331–1341, 2020\.
  * [10] J. Khazaei and D. H. Nguyen, “Distributed consensus for output power regulation of dfigs with on-site energy storage,” IEEE Transactions on Energy Conversion, vol. 34, no. 2, pp. 1043–1051, 2018.
  * [11] T. Meng, Z. Lin, and Y. A. Shamash, “Distributed cooperative control of battery energy storage systems in dc microgrids,” IEEE/CAA Journal of Automatica Sinica, vol. 8, no. 3, pp. 606–616, 2021.
  * [12] J. Cao, N. Schofield, and A. Emadi, “Battery balancing methods: A comprehensive review,” in 2008 IEEE Vehicle Power and Propulsion Conference, pp. 1–6, IEEE, 2008.
  * [13] Y. Zeng, Q. Zhang, Y. Liu, X. Zhuang, and H. Guo, “Hierarchical cooperative control strategy for battery storage system in islanded dc microgrid,” IEEE Transactions on Power Systems, 2021.
  * [14] M. Yazdanian and A. Mehrizi-Sani, “Distributed control techniques in microgrids,” IEEE Transactions on Smart Grid, vol. 5, no. 6, pp. 2901–2909, 2014.
  * [15] N. M. L. Tan, T. Abe, and H. Akagi, “Design and performance of a bidirectional isolated dc–dc converter for a battery energy storage system,” IEEE Transactions on Power Electronics, vol. 27, no. 3, pp. 1237–1248, 2011.
  * [16] G. Hu, “Robust consensus tracking of a class of second-order multi-agent dynamic systems,” Systems & Control Letters, vol. 61, no. 1, pp. 134–142, 2012.
  * [17] S. Shao, X. Liu, and J. Cao, “Prespecified-time synchronization of switched coupled neural networks via smooth controllers,” Neural Networks, vol. 133, pp. 32–39, 2021.


# Small scale quasi-geostrophic convective turbulence at large Rayleigh number

Tobias G. Oliver  Department of Physics, University of Colorado, Boulder, CO
80309, USA  Adrienne S. Jacobi  Department of Physics, University of Colorado,
Boulder, CO 80309, USA  Keith Julien  Department of Applied Mathematics,
University of Colorado, Boulder, CO 80309, USA  Michael A. Calkins  Department
of Physics, University of Colorado, Boulder, CO 80309, USA

###### Abstract

A numerical investigation of quasi-geostrophic Rayleigh-Bénard convection is
conducted in which the depth-invariant flows, including large-scale vortices,
are numerically suppressed. As the Rayleigh number $Ra$ is increased, the
Reynolds number and Nusselt number show evidence of approaching the diffusion-
free scalings of $Re\sim RaE/Pr$ and $Nu\sim Pr^{-1/2}Ra^{3/2}E^{2}$,
respectively, where $E$ is the Ekman number and $Pr$ is the Prandtl number.
For large $Ra$, the presence of depth-invariant flows, such as large-scale
vortices, yield heat and momentum transport scaling that exceeds those of the
diffusion-free scaling laws. The Taylor microscale does not vary significantly
with increasing $Ra$, whereas the integral length scale grows weakly. The
computed length scales remain $O(1)$ with respect to the linearly unstable
critical wavenumber and thus remain viscously controlled. We do not find a
point-wise Coriolis-Inertia-Archimedean (CIA) force balance in the turbulent
regime; interior dynamics are instead dominated by horizontal advection
(inertia), vortex stretching (Coriolis) and the vertical pressure gradient. A
secondary, sub-dominant balance between the buoyancy force and the viscous
force occurs in the interior and the ratio of the rms of these two forces is
found to approach unity with increasing $Ra$. This secondary balance is
attributed to the turbulent fluid interior acting as the dominant control on
the heat transport. These findings indicate that a pointwise CIA balance does
not occur in the high Rayleigh number regime of quasi-geostrophic Rayleigh-
Bénard convection. Instead, simulations are characterized by what may be
termed a non-local CIA balance in which the buoyancy force is dominant within
the thermal boundary layers and is spatially separated from the interior
Coriolis and inertial forces.

##  I Introduction

Convection plays an important role in the dynamics of many planets and stars,
where it serves as the power source for sustaining the magnetic fields of the
planets [29, 13, 24, 1] and the Sun [4]. Convection may also be a possible
driving mechanism for the observed large scale zonal winds [e.g. 12] and
vortices [e.g. 27] on the giant planets. The flows in these natural systems
are strongly forced and turbulent and can be constrained by the Coriolis
force. Studying rotationally constrained convective turbulence is therefore
important for improving understanding of such systems. However, experimental
[6, 5, 18, 32] and numerical [8, 11, 33] investigations have difficulty
accessing this parameter regime due to the extreme scale separation that
characterizes the dynamics. Asymptotic models play an important role in this
regard since they allow for significant computational savings by eliminating
physically unimportant dynamics while retaining the dominant force balance
that is thought to be representative of natural systems. In particular, the
asymptotic model for rapidly rotating convection in a planar geometry has been
used to advance our understanding of this system [16], and shows excellent
agreement with the results of DNS where an overlap of the parameter space is
possible [30, 22]. Nevertheless, the behavior of the system in the dual limit
of strong buoyancy forcing and rapid rotation is still not completely
understood [19]. In the present work we use numerical simulations of this
asymptotic model for investigating the scaling behavior of rotating convective
turbulence at previously inaccessible parameter values.

Rayleigh-Bénard convection, consisting of a fluid layer of depth $H$ confined
between plane parallel boundaries is a canonical system used for studying
buoyancy-driven flows. The two boundaries have temperature difference $\Delta
T$ and a constant gravitational field of magnitude $g$ points perpendicular to
the boundaries. The Rayleigh number

$Ra=\frac{g\alpha\Delta TH^{3}}{\nu\kappa},$ (1)

provides a non-dimensional measure of the buoyancy force. Here, $\alpha$ is
the thermal expansion coefficient, $\nu$ is the kinematic viscosity and
$\kappa$ is the thermal diffusivity. In most systems of interest, the flow is
strongly driven such that $Ra\gg 1$. In the presence of a background rotation
with angular frequency $\Omega$, the resulting convective dynamics are
considered rotationally constrained, or quasi-geostrophic (QG), when viscous
forces and inertia are both small relative to the Coriolis force. In non-
dimensional terms, quasi-geostrophy is characterized by small Rossby and Ekman
numbers, respectively defined by,

$Ro=\frac{U}{2\Omega H},\qquad E=\frac{\nu}{2\Omega H^{2}},$ (2)

where $U$ is a characteristic flow speed. In addition, the Reynolds number is
given by $Re=UH/\nu=Ro/E$; QG turbulence is characterized by $Re\gg 1$, and
therefore $E\ll Ro\ll 1$. As an example, for the Earth’s outer core estimates
suggest $Ro=O(10^{-7})$, $E=O(10^{-15})$ and $Re=O(10^{8})$ [e.g. 24].

An inverse kinetic energy cascade is known to occur in QG convection [17, 25,
7, 10, 19]. In the periodic plane layer geometry, the inverse cascade gives
rise to a depth-invariant large scale vortex (LSV) that tends to equilibrate
with lateral dimensions comparable to the horizontal size of the simulation
domain. The LSV can then dynamically influence the underlying small scale
convection. Previous studies have shown that the relative amount of kinetic
energy contained in the LSV and the small scale convection depends non-
monotonically on the forcing amplitude [10, 19]. In addition, the
characteristic speed of the LSV depends linearly on the aspect ratio of the
simulation domain [19, 21]. Although such coupling between the large and small
scale dynamics is directly relevant to natural systems, it nevertheless
complicates efforts to understand the strongly forced asymptotic regime of the
small scale convection.

Momentum and heat transport are charactrerised by $Re$ and the Nusselt number,
$Nu$, respectively. Understanding the scaling of these two quantities in the
high Rayleigh number regime of QG convection is vital for relating model
output to natural systems. QG convection dynamics are conveniently specified
by the asymptotic combination

$\widetilde{Ra}=RaE^{4/3},$ (3)

which we refer to as the reduced Rayleigh number. With this rescaling, the
onset of convection occurs at $\widetilde{Ra}\approx 8.7$, and the onset of
turbulence occurs when $\widetilde{Ra}\gtrsim 40$ [17]. Previous work has
found heat transport data that is consistent with the diffusion-free scaling,
$Nu\sim\widetilde{Ra}^{3/2}Pr^{-1/2}$, in the turbulent regime [15, 8], where
$Pr=\nu/\kappa$ is the Prandtl number. However, recent investigations that
allow for the amplitude of the LSV to fully saturate have found a stronger
scaling for $Nu$ [19], which might suggest that although the small scale
convection approaches a $Nu\sim\widetilde{Ra}^{3/2}Pr^{-1/2}$, the interaction
with the LSV provides an additional enhancement of heat transport. A similar
diffusion-free scaling for the flow speeds is given by $Re\sim RaE/Pr$
implying a rotational free-fall velocity $U_{r\\!f\\!f}=g\alpha\Delta
T/2\Omega$. Recent studies of convection in a spherical geometry support this
behavior [9]. However, similar to their heat transport findings, Ref. [19]
find that the Reynolds number scales more strongly than the diffusion-free
scaling.

A defining property of turbulence is the presence of a broad range of length
scales within the flow field. In spectral space this range of scales
translates into a broadband kinetic energy spectrum. The spectrum of non-
rotating turbulence broadens via a forward cascade in which there is a net
transfer of energy from low wavenumbers to high wavenumbers. A consequence of
this transfer of energy is that the length scale at which viscous dissipation
is dominant becomes ever smaller as $Re$ increases [e.g. 23].

The process by which spectral broadening occurs in QG convection is less
clear. In the limit $E\rightarrow 0$, linear theory shows that the onset of
convection occurs on a length scale of size $\ell=O\left(E^{1/3}\right)$. This
length scale arises because the viscous force facilitates convection by
simultaneously perturbing the geostrophic force balance and relaxing the
Taylor-Proudman constraint. Understanding which length scales emerge in the
strongly nonlinear regime is important for characterizing QG convective
turbulence. Previous studies suggest that the length scale varies with the
Rossby number as $\ell\sim Ro^{1/2}$ [e.g. 9], which is thought to arise from
the so-called Coriolis-Inertia-Archimedean (CIA) balance [14]. In terms of the
reduced Rayleigh number this scaling is equivalent to $\ell\sim
E^{1/3}\sqrt{\widetilde{Ra}/Pr}$ [e.g. 2]. Recent experimental work using
water as the working fluid has found a slightly weaker scaling than $\ell\sim
Ro^{1/2}$ [18]; this same investigation, and a numerical study of convection-
driven dynamos [33], finds that the correlation length scale of the vorticity
is approximately constant with increasing $\widetilde{Ra}$. It is presently
unknown whether this behavior persists in the limit of large $\widetilde{Ra}$.

In the present investigation we report on the results of numerical simulations
of the QG model of rotating convection in which the depth-invariant flows are
suppressed at each time step of the numerical integration. This suppression is
done to isolate the asymptotic behavior of the small-scale convection in the
absence of the LSV, and to simulate the largest accessible values of
$\widetilde{Ra}$ in an attempt to identify asymptotic trends in the scaling
behavior of various flow quantities. In Sec. II we provide an overview of the
QG model and numerical methods. Results and conclusions are given in Sec. III
and Sec. IV, respectively.

##  II Methods

###  II.1 Governing equations

In the present investigation we employ the asymptotically reduced form of the
governing equations given by Ref. [28]. When non-dimensionalized using the
small-scale viscous diffusion speed $\nu/\ell^{*}_{\nu}$, where
$\ell^{*}_{\nu}=HE^{1/3}$, these equations take the form

${\partial_{t}}\zeta+J[\ \psi,\zeta]\
-{\partial_{Z}}w=\nabla_{\perp}^{2}\zeta,$ (4) ${\partial_{t}}w+J[\ \psi,w]\
+{\partial_{Z}}\psi=\frac{\widetilde{Ra}}{Pr}\vartheta+\nabla_{\perp}^{2}w,$
(5) ${\partial_{t}}\vartheta+J[\ \psi,\vartheta]\
+w{\partial_{Z}}\overline{\Theta}=\frac{1}{Pr}\nabla_{\perp}^{2}\vartheta,$
(6)
${\partial_{Z}}\left(\overline{w\vartheta}\right)=\frac{1}{Pr}\partial_{Z}^{2}\overline{\Theta},$
(7)

where $t$ is time, the Cartesian coordinate system is denoted by $(x,y,Z)$,
the Jacobian operator is defined by $J[\ \psi,A]\
={\partial_{x}}\psi{\partial_{y}}A-{\partial_{y}}\psi{\partial_{x}}A$ for some
scalar field $A$, and the horizontal Laplacian operator is denoted by
$\nabla_{\perp}^{2}=\partial_{x}^{2}+\partial_{y}^{2}$. The vertical
components of vorticity and velocity are denoted by $\zeta$ and $w$,
respectively, $\psi$ is the geostrophic streamfunction, and $\vartheta$ is the
fluctuating temperature. The vorticity and streamfunction are related via
$\zeta=\nabla_{\perp}^{2}\psi$. The mean temperature is $\overline{\Theta}$,
where the overline denotes a horizontal average.

The governing equations are solved using a de-aliased pseudo-spectral method
in which the flow variables are expanded as Chebyshev polynomials in the
vertical dimension and Fourier series in the horizontal dimensions. A third-
order accurate implicit-explicit Runge-Kutta scheme is used to advance the
equations in time. Further details on the code can be found in Ref. [20]. The
geostrophic streamfunction is decomposed into depth-averaged and fluctuating
components (i.e. barotropic and baroclinic components, respectively); to
facilitate the removal of the LSV, the depth-averaged component is set to zero
at each time step of the numerical integration. Simulations are performed for
the range $10\leq\widetilde{Ra}\leq 280$ with $Pr=1$. Details of the numerical
simulations are provided in Table 1.

###  II.2 Diagnostic quantities

Heat and momentum transport are quantified with the Nusselt number $Nu$ and
the asymptotically reduced Reynolds number $\widetilde{Re}$, respectively. The
Nusselt number is defined by

$Nu=1+Pr\left<w\vartheta\right>,$ (8)

where the angled brackets denote a volume and time average. The asymptotically
rescaled Reynolds number is defined based on the vertical component of the
velocity field such that

$\widetilde{Re}=E^{1/3}Re=\sqrt{\left<w^{2}\right>}.$ (9)

For future reference, the relationship

$\frac{\widetilde{Ra}}{Pr^{2}}\left(Nu-1\right)=\varepsilon_{u},\quad
Nu=\varepsilon_{\theta}$ (10)

can be derived from the governing equations [17], where the viscous and
thermal dissipation rates are defined by, respectively,

$\varepsilon_{u}=\langle|\nabla_{\perp}w|^{2}\rangle+\langle\zeta^{2}\rangle,\quad\varepsilon_{\theta}=\langle|\nabla_{\perp}\theta|^{2}\rangle+\langle\left(\partial_{Z}\overline{\Theta}\right)^{2}\rangle.$
(11)

The integral scale, Taylor microscale, Kolmogorov scale, and temperature
integral scale are calculated for each simulation and defined by,
respectively,

$\ell_{I}=\frac{\int
k^{-1}\widehat{E}_{u}(k)dk}{\int\widehat{E}_{u}(k)dk},\quad\ell_{T}=\left(\frac{\int\widehat{E}_{u}(k)dk}{\int
k^{2}\widehat{E}_{u}(k)dk}\right)^{1/2},\quad\ell_{K}=\varepsilon_{u}^{-1/4},\quad\ell_{I}^{\vartheta}=\frac{\int
k^{-1}\widehat{E}_{\vartheta}(k)dk}{\int\widehat{E}_{\vartheta}(k)dk}$ (12)

where the time- and depth-averaged kinetic energy (temperature variance)
spectrum is denoted by $\widehat{E}_{u}(k)$
($\widehat{E}_{\vartheta}\left(k\right)$), and the horizontal wavenumber
vector is $\mathbf{k}=(k_{x},k_{y})$ with modulus
$k=\sqrt{k_{x}^{2}+k_{y}^{2}}$.

The Kolmogorov scale is often interpreted as the scale at which viscosity
dominates and the turbulent kinetic energy is dissipated into heat. The
integral scale and Taylor microscale are respectively interpreted as measures
of the correlation length of turbulent motions and the intermediate length
scale at which fluid viscosity significantly affects the dynamics of turbulent
motions. With these interpretations in mind, the relative ordering
$\ell_{I}>\ell_{T}>\ell_{K}$ holds. Furthermore, for rotating convection,
viscosity is a required ingredient in destabilizing a rotating fluid subject
to an adverse temperature gradient [3]. This implies that convective motions
are inherently influenced by viscosity and, as a consequence, the Taylor
microscale (assuming the physical interpetation holds) is tied to the linear
instability scale. i.e. we expect $\ell_{T}\sim 1$. The integral scale for the
temperature is interpreted as the scale at which the fluid is forced by
buoyancy.

$\widetilde{Ra}$ | $N_{x}\times N_{y}\times N_{Z}$ | $\Delta t$ | $n$ | $\widetilde{Re}$ | $Nu$ | $\ell_{I}$ | $\ell_{T}$ | $\ell_{K}$
---|---|---|---|---|---|---|---|---
$10$ | $384\times 384\times 96$ | $5\times 10^{-4}$ | $10$ | $0.679$ | $1.282$ | $1.074$ | $1.001$ | $0.917$
$12$ | $384\times 384\times 96$ | $5\times 10^{-4}$ | $10$ | $1.283$ | $1.791$ | $1.249$ | $1.029$ | $0.701$
$15$ | $384\times 384\times 96$ | $1\times 10^{-3}$ | $10$ | $2.335$ | $2.884$ | $1.628$ | $1.114$ | $0.555$
$20$ | $384\times 384\times 96$ | $5\times 10^{-4}$ | $10$ | $4.03$ | $5.215$ | $1.647$ | $1.101$ | $0.438$
$25$ | $384\times 384\times 96$ | $5\times 10^{-4}$ | $10$ | $5.776$ | $7.578$ | $1.643$ | $1.089$ | $0.371$
$30$ | $384\times 384\times 96$ | $5\times 10^{-4}$ | $10$ | $7.453$ | $10.02$ | $1.764$ | $1.094$ | $0.330$
$40$ | $384\times 384\times 96$ | $5\times 10^{-4}$ | $10$ | $10.19$ | $14.37$ | $1.753$ | $1.046$ | $0.283$
$50$ | $384\times 384\times 96$ | $5\times 10^{-4}$ | $10$ | $12.71$ | $18.54$ | $1.758$ | $1.005$ | $0.253$
$60$ | $384\times 384\times 96$ | $5\times 10^{-4}$ | $10$ | $14.79$ | $22.04$ | $1.739$ | $0.969$ | $0.231$
$70$ | $432\times 432\times 180$ | $1\times 10^{-4}$ | $10$ | $17.09$ | $26.25$ | $1.787$ | $0.956$ | $0.213$
$80$ | $576\times 576\times 192$ | $1\times 10^{-4}$ | $10$ | $19.43$ | $30.6$ | $1.866$ | $0.946$ | $0.199$
$100$ | $648\times 648\times 288$ | $2\times 10^{-4}$ | $10$ | $23.64$ | $38.79$ | $1.934$ | $0.929$ | $0.178$
$120$ | $648\times 648\times 324$ | $1\times 10^{-4}$ | $10$ | $28.41$ | $48.48$ | $2.106$ | $0.937$ | $0.161$
$140$ | $768\times 768\times 324$ | $1\times 10^{-4}$ | $10$ | $33.3$ | $60.22$ | $2.13$ | $0.933$ | $0.148$
$160$ | $720\times 720\times 360$ | $1\times 10^{-4}$ | $10$ | $38.14$ | $70.31$ | $2.235$ | $0.938$ | $0.137$
$180$ | $810\times 810\times 384$ | $5\times 10^{-5}$ | $10$ | $43.87$ | $85.87$ | $2.179$ | $0.916$ | $0.128$
$200$ | $960\times 960\times 450$ | $5\times 10^{-5}$ | $10$ | $47.9$ | $97.96$ | $2.261$ | $0.922$ | $0.121$
$220$ | $960\times 960\times 480$ | $5\times 10^{-5}$ | $10$ | $53.99$ | $114.0$ | $2.233$ | $0.904$ | $0.113$
$240$ | $1125\times 1125\times 540$ | $2\times 10^{-5}$ | $10$ | $59.65$ | $127.3$ | $2.582$ | $0.976$ | $0.108$
$260$ | $1200\times 1200\times 600$ | $2\times 10^{-5}$ | $10$ | $64.03$ | $140.9$ | $2.563$ | $0.962$ | $0.104$
$280$ | $1200\times 1200\times 675$ | $2\times 10^{-5}$ | $10$ | $70.88$ | $163.2$ | $2.551$ | $0.953$ | $0.098$
$100$ | $324\times 324\times 288$ | $2\times 10^{-4}$ | $5$ | $23.02$ | $37.47$ | $1.57$ | $0.898$ | $0.179$
$100$ | $480\times 480\times 288$ | $2\times 10^{-4}$ | $7$ | $23.83$ | $39.69$ | $1.774$ | $0.919$ | $0.177$
$100$ | $768\times 768\times 288$ | $2\times 10^{-4}$ | $12$ | $23.73$ | $39.30$ | $1.953$ | $0.925$ | $0.178$
$100$ | $972\times 972\times 288$ | $2\times 10^{-4}$ | $15$ | $23.82$ | $39.41$ | $1.984$ | $0.923$ | $0.177$
$100$ | $1296\times 1296\times 288$ | $2\times 10^{-4}$ | $20$ | $23.85$ | $39.51$ | $2.054$ | $0.932$ | $0.177$
$40$ | $192\times 192\times 96$ | $5\times 10^{-4}$ | $5$ | $10.18$ | $14.27$ | $1.548$ | $1.035$ | $0.285$
$40$ | $288\times 288\times 96$ | $5\times 10^{-4}$ | $7$ | $10.38$ | $14.69$ | $1.675$ | $1.048$ | $0.282$
$40$ | $480\times 480\times 96$ | $5\times 10^{-4}$ | $12$ | $10.29$ | $14.5$ | $1.886$ | $1.064$ | $0.283$
$40$ | $576\times 576\times 96$ | $5\times 10^{-4}$ | $15$ | $10.36$ | $14.63$ | $1.886$ | $1.06$ | $0.282$
$40$ | $768\times 768\times 96$ | $5\times 10^{-4}$ | $20$ | $10.32$ | $14.57$ | $1.86$ | $1.05$ | $0.282$
Table 1: Details of the numerical simulations in which depth-invariant flows are suppressed. The number of physical space grid points in each respective direction is specified by $N_{x}\times N_{y}\times N_{Z}$, the time step size is denoted by $\Delta t$, and the horizontal domain size, as specified by the integer number of critical horizontal wavelengths, is denoted by $n$. The asymptotically rescaled Reynolds number is $\widetilde{Re}$, $Nu$ is the Nusselt number, and the integral length scale, Taylor microscale and Kolmogorov scale are denoted by $\ell_{I}$, $\ell_{T}$ and $\ell_{K}$, respectively. All length scales are normalized by the critical wavelength $\lambda_{c}\approx 4.815$.  $\widetilde{Ra}$ | $N_{x}\times N_{y}\times N_{Z}$ | $\Delta t$ | $n$ | $\widetilde{Re}$ | $Nu$ | $\ell_{I}$ | $\ell_{T}$ | $\ell_{K}$
---|---|---|---|---|---|---|---|---
$20$ | $384\times 384\times 96$ | $5\times 10^{-4}$ | $10$ | $3.541$ | $4.01$ | $2.048$ | $1.189$ | $0.414$
$30^{*}$ | $128\times 128\times 64$ | $5\times 10^{-4}$ | $10$ | $7.219$ | $7.960$ | $--$ | $--$ | $--$
$40$ | $384\times 384\times 96$ | $5\times 10^{-4}$ | $10$ | $10.59$ | $11.79$ | $1.705$ | $1.023$ | $0.272$
$60^{*}$ | $256\times 256\times 96$ | $1\times 10^{-4}$ | $10$ | $16.82$ | $19.96$ | $--$ | $--$ | $--$
$80$ | $576\times 576\times 192$ | $1\times 10^{-4}$ | $10$ | $24.68$ | $30.92$ | $2.291$ | $0.928$ | $0.186$
$120$ | $648\times 648\times 324$ | $1\times 10^{-4}$ | $10$ | $41.4$ | $58.2$ | $2.805$ | $0.974$ | $0.140$
$160$ | $768\times 768\times 384$ | $5\times 10^{-5}$ | $10$ | $59.40$ | $98.06$ | $2.792$ | $0.919$ | $0.118$
$200$ | $960\times 960\times 450$ | $5\times 10^{-5}$ | $10$ | $84.21$ | $146.2$ | $3.032$ | $0.943$ | $0.098$
Table 2: Details of the numerical simulations for cases in which depth-
invariant flows are not suppressed. See Table 1 for specifics. Runs denoted
with an asterisk ($*$) are from Ref. [19].

##  III Results

###  III.1 Influence of simulation domain size

The presence of an inverse cascade complicates our understanding of rotating
convection since the flow speeds associated with the LSV are known to grow
linearly with the horizontal dimension of the simulation domain size [19].
This linear dependence is tied to the fact that the inverse cascade is halted
solely by the viscous force acting on the domain scale [21]. While we
eliminate the LSV (and all depth-invariant motions) in the current
investigation, it is nevertheless important to determine the domain size that
allows for convergence of key statistical quantities. Towards this end, a
series of simulations with $\widetilde{Ra}=40$ and $\widetilde{Ra}=100$ were
performed for varying domain sizes. Here we scale the horizontal dimension of
the simulation domain in integer multiples ($n$) of the critical horizontal
wavelength, $\lambda_{c}\approx 1.3048$. The effective resolution (i.e. number
of grid points per critical wavelength) was held approximately constant as the
domain size was varied. Figure 1 shows the sensitivity of $Nu$,
$\widetilde{Re}$, and the length scales $\ell_{I}$ and $\ell_{T}$ as a
function of $n$. We find that $Nu$, $\widetilde{Re}$, and $\ell_{T}$ show no
significant sensitivity to domain size with $n\geq 10$, which is consistent
with Ref. [19]. We find that the integral length scale $\ell_{I}$ exhibits
more sensitivity to the domain size. For $\widetilde{Ra}=40$, it appears that
$\ell_{I}$ converges to a nearly constant value for $n\geq 12$. For
$\widetilde{Ra}=100$ it is less clear whether convergence in $\ell_{I}$ is
achieved for the simulated range of $n$. The focus of the present
investigation is to reach the largest computationally affordable values of
$\widetilde{Ra}$; for this reason we choose $n=10$ for all of the data that is
presented in later subsections. Moreover, we find that the differences in
$\ell_{I}$ between simulations with $n=10$ and $n>10$ are less than $10\%$
over our investigated range of parameters.

(a)

(b)

(c)

Figure 1: Various flow quantities as a function of horizontal domain size (as
characterized by the number of critical horizontal wavelengths, $n$) for
$\widetilde{Ra}=100$ (◼) and $\widetilde{Ra}=40$ (▲). Error bars denote the
standard deviation of globally averaged quantities. LABEL:sub@sub:box_nondim
Nusselt number; LABEL:sub@sub:box_reynolds Reynolds number; (c) integral
length scale ($\ell_{I}$) and Taylor microscale ($\ell_{T}$).

###  III.2 Heat and momentum transport

(a)

(b)

Figure 2: Heat and momentum transport for all cases: (a) compensated Nusselt
number, $Nu/\widetilde{Ra}^{3/2}$; (b) compensated reduced Reynolds number,
$\widetilde{Re}/\widetilde{Ra}$. Blue data points are from simulations with a
non-zero barotropic component (LSV). Data with open symbols in panels (a) and
(b) shows the compensation required to flatten $Nu$ and $\widetilde{Re}$ for
the LSV data.

Panels (a) and (b) in Figure 2 show the compensated values
$Nu/\widetilde{Ra}^{3/2}$ and $\widetilde{Re}/\widetilde{Ra}$, respectively.
Both compensated quantities show steep initial increases up to
$\widetilde{Ra}\approx 30$, followed by a period of decrease before flattening
out with increasing $\widetilde{Ra}$. The cases with the LSV also show a
transition in scaling behavior near $\widetilde{Ra}\approx 30$, with an
additional transition in the heat transport at $\widetilde{Ra}\approx 60$ in
which a departure from the diffusion-free scaling is observed. At
$\widetilde{Ra}\approx 80$ we find a crossover in the heat transport; for
values less (more) than this the simulations without the LSV are more (less)
efficient at heat transport. The value of $\widetilde{Ra}\approx 80$ coincides
with the point at which the amplitude of the LSV reaches a maximum relative to
the small scale convection [19]. We also observe a crossover in the momentum
transport near $\widetilde{Ra}\approx 40$. The data suggests that in the
absence of the LSV, the baroclinic convection approaches scalings of
$Nu\sim\widetilde{Ra}^{3/2}$ and $\widetilde{Re}\sim\widetilde{Ra}$ for
$\widetilde{Ra}\gg 1$. For non-unity Prandtl numbers, these scalings extend to
the form $Nu\sim Pr^{-1/2}\widetilde{Ra}^{3/2}$ and
$\widetilde{Re}\sim\widetilde{Ra}/Pr$. However, we find that in comparison to
simulations without a LSV, simulations with a LSV show more efficient heat and
momentum transport for sufficiently large values of $\widetilde{Ra}$.

With the asymptotic scalings of $Nu\sim\widetilde{Ra}^{3/2}$ and
$\widetilde{Re}\sim\widetilde{Ra}$, equations (8) and (9) indicate that the
fluctuating temperature scales as $\vartheta\sim\widetilde{Ra}^{1/2}$. To test
this scaling we compute rms values of the fluctuating temperature over various
depths of the flow domain: (1) the maximum value taken at the ‘edge’, or
inflection point, of the thermal boundary layer; (2) averaged over the entire
layer depth $0\leq Z\leq 1$; and (3) averaged over the mid-depth range,
$0.4<Z<0.6$. Figure 3(a) shows these values; Figure 3(b) shows the compensated
values for the thermal boundary layer and entire depth. We find that both the
mid-depth and total depth rms values scale with $\widetilde{Ra}$ in the same
manner, whereas the thermal boundary layer shows a stronger dependence on
$\widetilde{Ra}$ (Figure 3LABEL:sub@temp_max_da), suggesting that the relative
thinness of the thermal boundary layer translates to an overall weaker
influence of this region on the observed scaling of the Nusselt number. From
the data, it was difficult to determine the powerlaw scaling for
$\vartheta_{rms}$ in the boundary layer, as results differed significantly
depending on the range of data used. For example, fitting the data in Figure
3LABEL:sub@temp_max_da for $\widetilde{Ra}>80$ gives a $\widetilde{Ra}^{0.78}$
relation while fitting over only the last three points ($\widetilde{Ra}>200$)
yields a $\widetilde{Ra}^{0.84}$ scaling. Nevertheless, we conclude that our
results for $\vartheta_{rms}$ in the thermal boundary layer are consistent
with the theoretically deduced $\widetilde{Ra}^{7/8}$ scaling discussed in
[15]. As previously reported in [28], we observe saturation of the interior
mean temperature gradient with a value of
$-\partial_{Z}\overline{\Theta}\approx 0.4$-$0.5$ (Figure 3LABEL:sub@dz_temp).

(a)

(b)

(c)

Figure 3: LABEL:sub@temp_da Various measures of the rms fluctuating
temperature vs. $\widetilde{Ra}$; LABEL:sub@temp_max_da compensated
fluctuating temperature; LABEL:sub@dz_temp mean temperature gradient vs.
$\widetilde{Ra}$. The vertical averaging range is indicated with subscripts in
(a) and (c).

###  III.3 Length scales and flow structure

(a)

(b)

(c)

(d)

Figure 4: Spectra for select values of $\widetilde{Ra}$: (a) kinetic energy
spectra vs. box normalized wavenumber; (b) compensated kinetic energy spectra;
(c) temperature variance spectra; (d) compensated temperature variance
spectra. Red dashed curves in show spectra for $\widetilde{Ra}=200$ with the
LSV.

Kinetic energy spectra, $\widehat{E}_{u}(k)$, and temperature variance
spectra, $\widehat{E}_{\vartheta}(k)$, are shown in Figure 4 for a range of
Rayleigh numbers. As expected, near the onset of convection for
$\widetilde{Ra}=10$, both $\widehat{E}_{u}(k)$ and
$\widehat{E}_{\vartheta}(k)$ show a well-defined peak near $k/k_{box}=10$,
corresponding to ten unstable wavelengths in the domain. A broadening of
dynamically active wavenumbers is observed in both spectra as $\widetilde{Ra}$
is increased. Figure 4LABEL:sub@energy_spectra_scaled shows the kinetic energy
spectra compensated by the Kolmogorov scaling of $k^{-5/3}$. We find that the
scaling of the kinetic energy spectra near $k/k_{box}\approx 10$ is somewhat
steeper than this classical scaling; a scaling closer to $k^{-5/2}$ provided a
better fit. This data suggests that an inertial subrange, at least in regards
to its classical definition of being absent of forcing and dissipation, does
not arise for the range of parameter values accessible here. The red dashed
curve shows the baroclinic kinetic energy spectra for a simulation with
$\widetilde{Ra}=200$ that includes the LSV. Comparison of the spectra at
$\widetilde{Ra}=200$ for cases with and without the LSV indicates that
significantly more energy is present in the smallest wavenumbers when the LSV
is present. This difference, along with the findings reported in Section
III.2, suggests that the presence of more energy in the largest length scales
contributes to more efficient heat and momentum transport. The temperature
variance spectra shown in Figure 4LABEL:sub@temp_spectra corroborate this
behavior; with the LSV we find significantly more energy in the $k/k_{box}=1$
mode in comparison to the same case without a LSV. We also note that in
comparison to $\widehat{E}_{u}(k)$, we find that $\widehat{E}_{\vartheta}(k)$
shows a significantly slower decay in amplitude with increasing $k$,
indicating that buoyancy forcing in the vicinity of the critical wavenumber
remains important even at very large values of $\widetilde{Ra}$. The
compensated temperature variance spectra shown in Figure 4(d) illustrates this
slow decay where an empirical scaling of $k^{-2/5}$ flattens the spectra in
the intermediate wavenumber range.

(a)

(b)

Figure 5: (a) Length scales versus $\widetilde{Ra}$ for all simulations where
$\ell_{I}$ is the integral scale, $\ell_{T}$ is the Taylor microscale,
$\ell_{K}$ is the Kolmogorov scale, and $\ell_{I}^{\vartheta}$ is the integral
scale for the fluctuating temperature. In (a) open symbols correspond to LSV
cases and the light blue line shows the $\widetilde{Ra}^{-1/2}$ scaling for
reference. (b) Comparison of computed length scales with linear stability
theory. $\widetilde{Ra}_{marg}$ is the marginal stability curve.
$\widetilde{Ra}_{fast}$ is the fastest growing mode in the linear stability
analysis about the saturated mean temperature profile. For large
$\widetilde{Ra}$, $\widetilde{Ra}_{fast}\sim k^{-8}$. The critical wavelength
and critical wavenumber are denoted by $(\ell_{c},k_{c})=(4.8154,1.3048)$.

Various length scales of the velocity field can be computed from the kinetic
energy spectra. Definitions for the integral length scale, the Taylor
microscale and the Kolmogorov length scale in terms of the kinetic energy
spectrum and vorticity (dissipation) are given in equation (12). Each of these
quantities are plotted in Figure 5LABEL:sub@lengthscales for all simulations.
We also plot the integral scale $\ell^{\vartheta}_{I}$ computed from
$\widehat{E}_{\vartheta}$, which represents a measure of the buoyancy driving
scale. For the LSV cases, the length scales are calculated using only the
baroclinic kinetic energy for better comparison with the present results. The
length scales are shown in units of the critical horizontal wavelength and we
find that all length scales are close to this critical wavelength just above
the onset of convection (i.e. when $\widetilde{Ra}=10$). Generally, the
integral length scale is found to be a slowly increasing function of
$\widetilde{Ra}$ regardless of whether the LSV is present or not. In
comparison to cases without the LSV, the LSV cases show a slightly faster
increase of $\ell_{I}$ for $\widetilde{Ra}>40$. At $\widetilde{Ra}=280$ we
find that the integral scale is approximately twice the size of the critical
wavelength. In comparison to the integral scale, the Taylor microscale shows
an initial increase up to $\widetilde{Ra}=15$, followed by a very slow
decrease. For $\widetilde{Ra}\geq 80,$ the average value of the Taylor
microscale is $\ell_{T}=0.94$ and the standard deviation is $0.02$. Both the
integral length scale and the Taylor microscale exhibit scaling behavior that
is weaker than the $\widetilde{Ra}^{1/2}$ scaling law. The Kolmogorov scale
exhibits an obvious decrease with $\widetilde{Ra}$ and, for the parameter
regime accessible here, is the only computed length scale that becomes
significantly different in value than the critical wavelength. For LSV cases,
$\ell_{K}$ is found to be slightly smaller in comparison to cases without the
LSV; this result is in agreement with the heat transport data which shows LSV
cases have larger heat transport and therefore larger viscous dissipation. The
integral temperature scale remains $O\left(1\right)$ up to
$\widetilde{Ra}=280$, which suggests that the buoyancy forcing scale occurs at
a viscous length scale. Around $\widetilde{Ra}\geq 200$,
$\ell_{I}^{\vartheta}$ begins to slowly increase. This behaviour is co-
incident with a decrease in $\ell_{I}^{\vartheta}$ for the LSV cases. It
appears that a crossover in $\ell_{I}^{\vartheta}$ for cases with the LSV and
cases without a LSV occurs around $\widetilde{Ra}=200$. This result suggests
that at large $\widetilde{Ra}$, the LSV facilitates the buoyancy forcing at
smaller scales than those without the LSV.

Due to the similarity between the computed length scales and the linearly
unstable wavelength, we show an alternative presentation of length scales
based on comparisons with the marginal stability boundary and the fastest
growing linearly unstable modes in Figure 5LABEL:sub@fastest_mode. The
marginal stability boundary is defined by

$\widetilde{Ra}_{marg}=k^{4}+\pi^{2}k^{-2}.$ (13)

The small and large wavenumber scalings of $k^{-2}$ and $k^{4}$ are also shown
where we note that the former is consistent with CIA theory in that it
represents a diffusion-free scaling behavior. We find that none of the
computed length scales exhibit scaling behavior that is similar to this
diffusion-free trend.

For any constant (non-zero) mean temperature gradient, a linear stability
analysis can be performed on equations (4)-(7) to acquire a growth rate
$\lambda\left(k,\widetilde{Ra}\right).$ The fastest growing mode is found by
taking the global maximum of $\lambda$ for a given value of $\widetilde{Ra}.$
In the limit of large $\widetilde{Ra}$, this analysis suggests that the length
scale associated with the fastest growing mode scales like
$\ell_{fast}\sim\widetilde{Ra}^{-1/8}$; we find that $\ell_{I}$ shows a
scaling behavior that is similar to this trend.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6:  Snapshots of the fluctuating temperature: (a,d)
$\widetilde{Ra}=40$; (b,e) $\widetilde{Ra}=160$; (c,f) $\widetilde{Ra}=280$.
Integral temperature scale is given by the black bars. (a)-(c) Top row:
$x$-$Z$ slices; bottom row: $x$-$y$ midplane slices. The units of $x$ and $y$
are given in number of critical wavelengths. Dot-dashed lines mark the
intersection of the two planes ($x$-$Z$ and $x$-$y$). (d)-(f) Horizontal cross
sections within the (upper) thermal boundary layer.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7:  Snapshots of the vorticity: (a,d) $\widetilde{Ra}=40$; (b,e)
$\widetilde{Ra}=160$; (c,f) $\widetilde{Ra}=280$. Integral scale, Taylor
microscale and ($10\times$) Kolmogorov scale are given by the black bars.
Green corresponds to cyclonic and purple to anti-cyclonic.

Instantaneous physical space visualizations of the fluctuating temperature and
vertical vorticity are shown in Figures 6 and 7, respectively, for
$\widetilde{Ra}=40$ (a,d), $\widetilde{Ra}=160$ (b,e) and $\widetilde{Ra}=280$
(c,f). For each value of $\widetilde{Ra}$, all visualizations are taken at the
same instant in time. Vertical slices are shown in the top row of panels
(a)-(c) and horizontal slices taken from the mid-plane of the fluid layer are
shown in the bottom row of panels (a)-(c); the dot-dashed lines in (a)-(c)
mark the intersection of the vertical and horizontal planes in the top and
bottom rows, respectively. Panels (d)-(f) show horizontal slices at depths in
which $\vartheta_{rms}$ reaches its maximum value, i.e. within the thermal
boundary layer. The visualizations confirm that the characteristic size of the
large scale flow patterns, as quantified by the integral length scale, is
weakly dependent on the Rayleigh number. Moreover, these large scale
structures show significant axial coherence even in the absence of depth-
averaged flows. Even for $\widetilde{Ra}=280$, flow structures that span the
entire depth can be identified. However, it can be observed in both fields
that the vertical length scale of coherence within the interior decreases with
increasing $\widetilde{Ra}$. This is discussed in the context of non-local CIA
theory in the final section. The thermal boundary layer visualizations show
evidence of strong spatial correlation between the thermal and vortical
fields. While the size of the largest vortical structures remain largely
insensitive to $\widetilde{Ra}$, it is observed that the filamentary
structures become increasingly finer which is consistent with the theoretical
findings in Ref. [15]. For the vortical field, this fine scale is imprinted
into the interior whereas the interior spatial structure of the temperature
field is observed to be de-correlated with the boundary. In the next section,
we associate this effect to the thermal field behaving as an advective-
diffusive scalar stirred by the vortical field and thermally dissipated.

###  III.4 Balances

(a)

(b)

(c)

Figure 8: Logarithmic vertical profiles of rms values of each term in the
governing equations for $\widetilde{Ra}=200$: (a) vertical vorticity; (b)
vertical momentum; (c) fluctuating temperature.

In an effort to understand the scaling behavior of the various quantities
discussed thus far, we compute rms values of all terms appearing in equations
(4)-(6). Vertical profiles of these quantities are shown in Figure 8 for
$\widetilde{Ra}=200$, which is a value that is representative of the turbulent
regime. The data shown in each plot is computed by squaring each term in the
respective equations, averaging this quantity over the horizontal plane,
taking the square root, and then averaging in time. The material derivative
$D_{t}$ is the rms of the sum of all forcing terms with the advective terms
subtracted off (i.e. for the vorticity equation (4), rms($D_{t}\zeta$) =
rms($\partial_{Z}w+\nabla_{\perp}^{2}\zeta$)). Note that a logarithmic scale
is used on the vertical axis of each figure to illustrate the differences in
balances that occur in the interior with those that occur within the thermal
boundary layer. We find that advection and the time derivative terms are
largest in each of the three equations throughout the fluid layer, which
agrees with the results of Ref. [19] for simulations in which depth-invariant
flows were present. However, the dynamics are controlled by the material
derivative, rather than the time derivative and advection separately. Panel
(a) shows that $D_{t}\zeta$ and $\partial_{z}w$ are nearly identical in
magnitude throughout the fluid layer. However, while smaller, the diffusion of
vorticity remains comparable in magnitude to these two terms; in the fluid
interior ($0.5>Z>0.1$) we find rms($\nabla_{\perp}^{2}\zeta$) $\sim 2\times
10^{3}$, whereas rms($D_{t}\zeta$) $\sim 3\times 10^{3}$ and
rms($\partial_{z}w$) $\sim 3\times 10^{3}$. Within the thermal boundary layer
($Z\lesssim 10^{-3}$) we find that all terms in the vorticity equation are
important, indicating that the large amplitude vortical motions are directly
influenced by viscous diffusion at the boundaries.

For the vertical momentum equation balances shown in Figure 8(b) we find that
$D_{t}w$ and $\partial_{Z}\psi$ are of nearly identical magnitude in the fluid
interior, whereas the buoyancy force and diffusion terms are comparable in
magnitude and the smallest of all terms in the interior. These results show
that the vertical pressure gradient acts as the dominant driver of vertical
motion in the interior. In the thermal boundary layer we find that all terms
in the vertical momentum equation become comparable in magnitude, though as a
consequence of impenetrability, diffusion remains the smallest of all terms.
Figure 8(c) shows a tendency in the interior for the horizontal material
advection of fluctuating temperature to be balanced by horizontal thermal
diffusion, whereas all terms become comparable in magnitude within the thermal
boundary layer.

Figure 9: Select ratios of rms forces from vertical momentum and vorticity
equations. In the vertical momentum equation, the ratio of diffusion to
buoyancy approaches unity. The ratio of diffusion to pressure gradient (vortex
stretching) in the vertical momentum (vorticity) equation is a slowly
decreasing function of $\widetilde{Ra}$ for large $\widetilde{Ra}$.

Figure 9 show ratios of rms values of various terms in the vertical momentum
and vertical vorticity equations as function of $\widetilde{Ra}$. The rms
values are the depth-averaged values of vertical profiles similar to those
shown in Figure 8(b). We find that the ratio of the vertical viscous force to
the buoyancy force approaches unity as $\widetilde{Ra}$ is increased. Given
the scaling behavior of the rms fluctuating temperature, $\vartheta_{rms}$ (as
shown in Figure 3(a)) this data indicates that the viscous force increases at
a rate faster than $\widetilde{Ra}^{1/2}$ until unity is reached in the force
balance. This balance between the viscous force and the buoyancy force is
attributed to the observation that the fluid interior controls heat transport
in the large Rayleigh number regime of rapidly rotating convection [15]. The
ratios of the viscous force to the pressure gradient force and the diffusion
of vertical vorticity to vortex stretching are slowly decreasing functions of
$\widetilde{Ra}$. From our data, while it appears that the magnitude of the
globally averaged influence of viscosity is subdominant from that of the
pressure gradient or vortex stretching in the limit of large $\widetilde{Ra}$,
it is unclear whether asymptotic subdominance holds for
$\widetilde{Ra}\rightarrow\infty$.

##  IV Discussion and Conclusions

The QG equations (4)–(7) represent the asymptotic low Rossby number limit of
the buoyantly forced Navier-Stokes equations. Simulations of these equations
were performed in which the depth-invariant flows were entirely suppressed at
each numerical time step. The suppression was done to isolate the dynamics of
the small scale convection and to enable simulations at previously
inaccessible values of $\widetilde{Ra}$. A comparison was made with data from
previously published simulations of the asymptotic model in which large
amplitude, depth-invariant flows (e.g. LSVs) were present. This comparison has
allowed for additional insight into the physics of small scale, QG convective
turbulence, which drives the inverse kinetic energy cascade in this system.
Asymptotically reduced Rayleigh and Reynolds numbers up to
$\widetilde{Ra}=280$ and $\widetilde{Re}\approx 70$ were simulated, which
represents the most extreme parameter regime accessed to date for QG
convection in either the plane layer or spherical geometry. Recent DNS studies
have reached small scale Reynolds numbers up to $\widetilde{Re}\approx 33$ in
the plane layer [34, 33], and up to $\widetilde{Re}\approx 24$ in a spherical
shell geometry [26]. Ref. [9] used a non-asymptotic two-dimensional model in
the equatorial region of a full sphere to reach $\widetilde{Re}\approx 20$.
For context, estimates suggest that characteristic flow speeds in the Earth’s
outer core yield reduced Reynolds numbers of $\widetilde{Re}\approx
10^{3}-10^{4}$, clearly indicating that it is necessary to understand the
physics of QG convection at large values of the reduced parameters.

A subset of simulations were performed in which the horizontal dimensions of
the simulation domain were systematically varied. In agreement with previous
studies of non-rotating convection [e.g. 31], we find that different
quantities show a different dependence on the simulation domain size (Figure
1). Globally averaged statistics, such as the Nusselt number and the Reynolds
number, converge rapidly for horizontal dimensions that are larger than ten
critical horizontal wavelengths. The Taylor microscale remains approximately
constant across the different domain sizes that were simulated. The integral
length scale of the flow field converges to a constant value with increasing
domain size, but at a rate that may be influenced by the particular value of
$\widetilde{Ra}$. However, the relative differences observed are small for our
investigated range of $\widetilde{Ra}$; a conclusion that is supported by the
rapid convergence of the Reynolds number with increasing domain size.

Our investigation finds evidence that in the absence of depth-invariant flows
both the heat and momentum transport approach asymptotic scaling regimes for
large Rayleigh numbers (Figure 2). When generalized to non-unity Prandtl
numbers, the data is well represented by scalings of the form $Nu\sim
Pr^{-1/2}\widetilde{Ra}^{3/2}$ and $\widetilde{Re}\sim\widetilde{Ra}/Pr$. When
recast in terms of large scale quantities these scalings become $Nu\sim
Pr^{-1/2}Ra^{3/2}E^{2}$ and $Re\sim RaE/Pr$, which represent diffusion-free
scalings when interpreted on the large domain scale [2]. The simulations that
include depth-invariant flows exhibit heat and momentum transport that is more
efficient than the aforementioned diffusion-free scalings [19]. This diffusion
dependence is of interest when considering applications to geophysical and
astrophysical systems in which large-scale flows, such as zonal jets and
vortices, are present.

Kinetic energy spectra and temperature variance spectra were computed for all
values of $\widetilde{Ra}$ (Figure 4). We do not find evidence of a trend
toward a Kolmogorov-like inertial subrange in which the kinetic energy spectra
scales with the horizontal wavenumber as $k^{-5/3}$; a stronger dependence
closer to $k^{-5/2}$ is observed. The kinetic energy spectra show that there
is a build up of energy at smaller wavenumbers as $\widetilde{Ra}$ is
increased, though the energy contained in these larger scale structures
remains comparable to that contained in scales similar to the critical
wavelength. The temperature variance spectra exhibit similar behavior,
indicating that the system remains forced on a length scale that is comparable
to the critical wavelength at the highest Rayleigh number investigated.

The evolution of various length scales in the system was investigated for
varying $\widetilde{Ra}$. We find that the integral length scale $\ell_{I}$
increases with $\widetilde{Ra}$, but at a rate that is slower than the
diffusion free scaling $\propto(\widetilde{Ra}/Pr)^{1/2}$ suggested by a CIA
balance. The Taylor microscale remains nearly constant with increasing
$\widetilde{Ra}$, whereas the Kolmogorov length scale decreases with
increasing $\widetilde{Ra}$. The temperature integral scale
$\ell_{I}^{\vartheta}$ was also found to remain approximately constant across
the range $20\lesssim\widetilde{Ra}\lesssim 200$, whereas a slow increase was
observed for $\widetilde{Ra}>200$. One of the key findings of this
investigation is that all of these length scales remain comparable to the
linearly unstable critical wavelength that emerges at the onset of convection.
We note that the scaling behavior of the integral scale and Taylor microscale
are consistent with recent laboratory experiments of rotating convection [18]
and a numerical study of rotating convection-driven dynamos in the plane layer
[33]. These findings provide evidence that the broadening of the length scales
in rotating convective turbulence is an extremely slow function of the
Rayleigh number, and occurs in a fundamentally different manner in comparison
to non-rotating turbulence. In particular, studies of non-rotating convection
find that both the integral length scale and the Taylor microscale decrease
rapidly with increasing Rayleigh number [e.g. 35].

An explanation for the aforementioned contrast to the non-rotating case
resides in the fact that viscosity is a required ingredient in destabilizing a
rotating fluid subject to an adverse temperature gradient in the presence of
the Taylor-Proudman constraint [3]. Convective motions are therefore
inextricably influenced by viscosity and, by definition, the Taylor microscale
is tied to the linear instability scale $\ell^{*}_{\nu}=E^{1/3}H$. This
connection to linear theory is demonstrated by comparisons between the
computed lengths and the marginal stability curve (Figure
5LABEL:sub@fastest_mode). It is found that the dependence of the integral
length scale with increasing $\widetilde{Ra}$ follows the fastest growing,
viscosity dependent, linear mode $\propto\widetilde{Ra}^{1/8}$, rather than
dissipation-free scale predicted by CIA theory $\propto\widetilde{Ra}^{1/2}$.
We note that the latter scaling resides at the long wavelength marginal
stability boundary for convective onset. The comparisons with linear theory in
which a mean temperature profile with gradient
$\partial_{Z}\overline{\Theta}=-1$ might remain pertinent in the nonlinear
regime given that the observed nonlinear profile saturates with a gradient
$\sim-0.5$ (Figure 3LABEL:sub@dz_temp). The permissible convectively unstable
length scales within the bulk thus adhere to the same analytic structure as
deduced from the marginal stability at linear onset.

Congruent with the lack of agreement between the measured integral length
scale and that predicted by CIA theory, the force balances also do not show
evidence of a CIA balance within the fluid interior, as might be predicted
from scaling arguments. Instead we find a trend indicating that this balance
is never approached as the system becomes more strongly forced since there is
a subdominant balance between the buoyancy force and viscous diffusion of
vertical momentum in the interior. The ratio of the rms of these two forces
approaches unity as $\widetilde{Ra}$ is increased, which we attribute to the
fact that the fluid interior controls heat transport in rapidly rotating
convection [e.g. 15]. The buoyancy force is strongest within the thermal
boundary layer and comparable to all other forces, including viscous forces,
in that region. Moreover, the strong, horizontal, vortical flows that are
present at the boundary are balanced entirely by viscosity. Convective
overturning motions in the interior are driven primarily by vertical pressure
gradients, and vortical motions are driven by vortex-stretching. However,
neither of these forces are sources of potential energy injection which must
then be provided from within the boundary regions.

The observed juxtaposition between global heat and momentum transport laws
that support a dissipation-free CIA balance and an observed integral length
scaling that does not can be resolved by analyzing force balances within the
fluid interior and boundary layers. Indeed, balances in the governing
equations can often be used to derive various scaling laws for the Rayleigh
number dependence of flow length scales. Towards this end we have computed the
magnitudes of the various terms appearing in the governing equations of the
asymptotic model (Figure 8). In agreement with previous simulations with
depth-invariant flows, we find that within the interior the turbulent regime
is predominantly characterized by the passive horizontal advection of each of
the flow quantities (vertical vorticity, vertical velocity and fluctuating
temperature). Details of the evolution of the flow quantities are best viewed
by following fluid elements in a Lagrangian framework. This captures the
hierarchy of secondary forces which drive the horizontal material advection,
as defined by $D_{t}$. We find that the following analysis is best pursued
utilizing a dimensional approach, followed by a reconversion to dimensionless
quantities. Considering the vertical momentum and vorticity equations in
dimensional form, we find

$\displaystyle D^{*}_{t^{*}}w^{*}\sim\partial_{Z^{*}}p^{*}$
$\displaystyle>g\alpha\theta^{*}\sim\nu\nabla^{*2}_{\perp}w^{*},$ (14)
$\displaystyle D^{*}_{t^{*}}\zeta^{*}\sim 2\Omega\partial_{Z^{*}}w^{*}$
$\displaystyle>\nu\nabla^{*2}_{\perp}\zeta^{*}.$ (15)

We therefore find an interior Coriolis-Inertial (CI) balance in which the
geostrophic pressure gradient and vortex stretching forces are the dominant
secondary forces in the hierarchy. Furthermore, the rms buoyancy and rms
viscous forces are comparable in magnitude and are the smallest of all terms.
We note that these balances are consistent with those found in DNS studies of
rotating convection [11] and rotating convection-driven dynamos [33]. They
imply that interior forcing and viscous dissipation play a subdominant role in
the dynamics. Consequently, the observed dynamics must be driven externally by
buoyant (A)rchimedean forces arising within the thermal boundary layer.

To summarize, we observe what may be termed a non-local CIA balance where the
action of the Archimedean buoyancy force occurring within the thermal boundary
layers is spatially separated from the interior Coriolis and inertial forces.
It can now be shown that this finding resolves the dichotomy of not finding a
diffusion-free integral scale $\ell^{*}_{I}$ as predicted by a local CIA
balance. Assuming velocities achieve rotational free-fall, as presently
observed, and that pressure is the geostrophic streamfunction, i.e.,

$\mathbf{u}^{*}\sim U_{r\\!f\\!f}=\frac{g\alpha\Delta T}{2\Omega},\quad
p^{*}\sim 2\Omega U_{r\\!f\\!f}\ell^{*}_{I},\quad\zeta^{*}\sim
U_{r\\!f\\!f}/\ell^{*}_{I},$ (16)

it follows from the interior CI balance given in (14,15) that

$\partial_{Z^{*}}\sim\frac{1}{h^{*}},\quad D^{*}_{t^{*}}\sim
2\Omega\frac{\ell^{*}_{I}}{h^{*}}.$ (17)

Here $h^{*}<H$ is the vertical correlation height of interior convective
motions. Viscous processes are estimated via
$\nabla^{*2}_{\perp}\sim(1/\ell^{*}_{d})^{2}$ where $\ell^{*}_{d}$ denotes the
(optimal) dissipation length scale. With these findings the sub-dominance and
equivalence of buoyancy and viscous forcing in (14,15) imply respectively

$h^{*}<\frac{2\Omega}{\nu}\ell^{*}_{I}\ell^{*2}_{d},\qquad\frac{\theta^{*}}{\Delta_{T}}\sim\frac{\nu}{2\Omega\ell^{*2}_{d}}$
(18)

We note that $\theta^{*}=E^{1/3}\Delta_{T}\theta$ in the QG regime and that
the boundary forced interior must also induce non-dimensional temperature
fluctuations $\theta\sim\left(\widetilde{Ra}/Pr\right)^{1/2}$ capable of
transporting the dissipation-free heat flux in a manner consistent with the
exact relation for the thermal dissipation rate (11), (see Figure
3LABEL:sub@temp_da). Equation (18) can be reformulated non-dimensionally as

$h<\ell_{I}\ell^{2}_{d},\quad\ell_{d}=\theta^{-1/2}\sim\left(\frac{\widetilde{Ra}}{Pr}\right)^{-1/4}$
(19)

with non-dimensional length scales $h=h^{*}/H$,
$\ell_{I}=\ell^{*}I/\ell^{*}_{\nu}$ and $\ell_{d}=\ell^{*}d/\ell^{*}_{\nu}$,
where $I$ and $d$ are non-dimensional coefficients; recall
$\ell^{*}_{\nu}=E^{1/3}H$. Most importantly, we note that the absence of a
local Archimedean force in the interior balance places no dissipation-free
restrictions on the scaling of the injection scale which is consistent with
the simulated observation $\ell_{I}\propto\widetilde{Ra}^{1/8}$. It now
follows that $h\leq O(\widetilde{Ra}^{-3/8})$ indicating a shortening of
vertical correlation length with increasing $\widetilde{Ra}$ as evident in
Figure 6 and 7. The interior force balances for the fluctuating temperature
obeys

$D^{\theta*}_{t^{*}}\theta^{*}\sim
w^{*}\partial_{Z^{*}}\overline{\Theta}^{*}\sim\kappa\nabla^{*2}_{\perp}\theta^{*}.$
(20)

such that unlike the momentum field, material advection (or stirring) of the
temperature field is balanced by dissipation (mixing). This can observed by
comparing the spatial morphologies of the midplane snapshots for temperature
and vorticity. The former exhibiting broader structures due to enhanced
mixing.

The results presented in this study highlight the non-trivial behavior of
rapidly rotating convection and provide a foundation for comparison with DNS
and experiment in both the plane-layer and spherical geometries. By removal of
the LSV, our study was able to investigate the origin of the inverse cascade
and our results demonstrate how rotating convection departs from the theory of
isotropic, homogenous turbulence. Most significantly, we have shown that
diffusion free force-balance arguments, even in the regime of large $Re$
provide an incomplete picture of the convective dynamics, and are insufficient
when attempting to identify the pertinent length scales in the flow. Instead,
our investigations reveal that the interior, despite controlling the global
momentum and heat transport, is forced externally by the boundary layers.
Therefore, we determine that QG dynamics at large Rayleigh number should not
be considered diffusion free.

Several open questions are apparent from this investigation, including: Does
the geostrophic regime accessed in this study accurately represent the flow
regime for $\widetilde{Ra}\to\infty$? If not, does a new, higher
$\widetilde{Ra}$ regime exist in which the impact of molecular dissipation is
diminished? To what extent does domain geometry impact the ultimate scaling
theory? Investigating these questions remains challenging to laboratory
studies and DNS given present difficulties in investigating broad ranges of
the extreme geostrophic parameter regime.

## Declaration of Interests

The authors report no conflict of interest.

## Acknowledgements

The authors gratefully acknowledge funding from the National Science
Foundation through grants EAR-1945270 (TGO, MAC), SPG-1743852 (TGO, MAC), and
DMS-2009319 (KJ). The simulations were conducted on the Summit and Stampede2
supercomputers. Summit is a joint effort of the University of Colorado Boulder
and Colorado State University, supported by NSF awards ACI-1532235 and
ACI-1532236. Stampede2 was made available through Extreme Science and
Engineering Discovery Environment (XSEDE) allocation PHY180013.

## References

  * Aurnou et al. [2015] Aurnou, J. M., Calkins, M. A., Cheng, J. S., Julien, K., King, E. M., Nieves, D., Soderlund, K. M. & Stellmach, S. 2015 Rotating convective turbulence in Earth and planetary cores. Phys. Earth Planet. Int. 246, 52–71.
  * Aurnou et al. [2020] Aurnou, J. M., Horn, S. & Julien, K. 2020 Connections between nonrotating, slowly rotating, and rapidly rotating turbulent convection transport scalings. Phys. Rev. Res.  2 (4), 043115\.
  * Chandrasekhar [1961] Chandrasekhar, S. 1961 Hydrodynamic and Hydromagnetic Stability. U.K.: Oxford University Press.
  * Charbonneau [2014] Charbonneau, P. 2014 Solar dynamo theory. Annu. Rev. Astron. Astrophys. 52, 251–290.
  * Cheng et al. [2015] Cheng, J. S., Stellmach, S., Ribeiro, A., Grannan, A., King, E. M. & Aurnou, J. M. 2015 Laboratory-numerical models of rapidly rotating convection in planetary cores. Geophys. J. Int. 201, 1–17.
  * Ecke & Niemela [2014] Ecke, R. E. & Niemela, J. J. 2014 Heat transport in the geostrophic regime of rotating Rayleigh-Bénard convection. Phys. Rev. Lett. 113 (114301).
  * Favier et al. [2014] Favier, B., Silvers, L. J. & Proctor, M. R. E. 2014 Inverse cascade and symmetry breaking in rapidly rotating Boussinesq convection. Phys. Fluids 26 (096605).
  * Gastine et al. [2016] Gastine, T., Wicht, J. & Aubert, J. 2016 Scaling regimes in spherical shell rotating convection. J. Fluid Mech. 808, 690–732.
  * Guervilly et al. [2019] Guervilly, C., Cardin, P. & Schaeffer, N. 2019 Turbulent convective length scale in planetary cores. Nature 570 (7761), 368–371.
  * Guervilly et al. [2014] Guervilly, C., Hughes, D. W. & Jones, C. A. 2014 Large-scale vortices in rapidly rotating Rayleigh-Bénard Convection. J. Fluid Mech. 758, 407–435.
  * Guzmàn et al. [2021] Guzmàn, A., M., Madonia, J., Cheng, R., Ostillo-Mònico, H., Clercx & R., Kunnen 2021 Force balance in rapidly rotating rayleigh-bènard convection. J. Fluid Mech. 928, A16.
  * Heimpel et al. [2016] Heimpel, M., Gastine, T. & Wicht, J. 2016 Simulation of deep-seated zonal jets and shallow vortices in gas giant atmospheres. Nature Geoscience 9 (1), 19.
  * Jones [2011] Jones, C. A. 2011 Planetary magnetic fields and fluid dynamos. Annu. Rev. Fluid Mech. 43, 583–614.
  * Jones [2015] Jones, C. A. 2015 Thermal and compositional convection in the outer core. In Treatise on Geophysics (ed. G. Schubert), , vol. 8, pp. 115–159. Elsevier.
  * Julien et al. [2012a] Julien, K., Knobloch, E., Rubio, A. M. & Vasil, G. M. 2012a Heat transport in Low-Rossby-number Rayleigh-Bénard Convection. Phys. Rev. Lett. 109, 254503.
  * Julien et al. [1998] Julien, K., Knobloch, E. & Werne, J. 1998 A new class of equations for rotationally constrained flows. Theoret. Comput. Fluid Dyn. 11, 251–261.
  * Julien et al. [2012b] Julien, K., Rubio, A. M., Grooms, I. & Knobloch, E. 2012b Statistical and physical balances in low Rossby number Rayleigh-Bénard convection. Geophys. Astrophys. Fluid Dyn. 106 (4-5), 392–428.
  * Madonia et al. [2021] Madonia, M., Guzmán, A. J. Aguirre, Clercx, H. J. H. & Kunnen, R. P. J. 2021 Velocimetry in rapidly rotating convection: Spatial correlations, flow structures and length scales (a). Europhys. Lett. 135 (5), 54002.
  * Maffei et al. [2021] Maffei, S., Krouss, M. J., Julien, K. & Calkins, M. A. 2021 On the inverse cascade and flow speed scaling behaviour in rapidly rotating Rayleigh–Bénard convection. J. Fluid Mech.  913, A18.
  * Marti et al. [2016] Marti, P., Calkins, M. A. & Julien, K. 2016 A computationally efficent spectral method for modeling core dynamics. Geochem. Geophys. Geosys. 17 (8), 3031–3053.
  * Nicoski et al. [2022] Nicoski, J. A., Yan, M. & Calkins, M. A. 2022 Quasistatic magnetoconvection with a tilted magnetic field. Phys. Rev. Fluids  7 (4), 043504.
  * Plumley et al. [2016] Plumley, M., Julien, K., Marti, P. & Stellmach, S. 2016 The effects of Ekman pumping on quasi-geostrophic Rayleigh-Bénard convection. J. Fluid Mech.  803, 51–71.
  * Pope [2000] Pope, S. B. 2000 Turbulent Flows. Cambridge: Cambridge University Press.
  * Roberts & King [2013] Roberts, P. H. & King, E. M. 2013 On the genesis of the Earth’s magnetism. Rep. Prog. Phys. 76 (096801).
  * Rubio et al. [2014] Rubio, A. M., Julien, K., Knobloch, E. & Weiss, J. B. 2014 Upscale energy transfer in three-dimensional rapidly rotating turbulent convection. Phys. Rev. Lett. 112, 144501\.
  * Schaeffer et al. [2017] Schaeffer, N., Jault, D., Nataf, H.-C. & Fournier, A. 2017 Geodynamo simulations with vigorous convection and low viscosity. Geophys. J. Int., https://arxiv.org/abs/1701.01299 , arXiv: 1701.01299.
  * Siegelman et al. [2022] Siegelman, L., Klein, P., Ingersoll, A. P., Ewald, S. P., Young, W. R., Bracco, A., Mura, A., Adriani, A., Grassi, D., Plainaki, C. & others 2022 Moist convection drives an upscale energy transfer at jovian high latitudes. Nature Physics 18 (3), 357–361.
  * Sprague et al. [2006] Sprague, M., Julien, K., Knobloch, E. & Werne, J. 2006 Numerical simulation of an asymptotically reduced system for rotationally constrained convection. J. Fluid Mech. 551, 141–174.
  * Stanley & Glatzmaier [2010] Stanley, S. & Glatzmaier, G. A. 2010 Dynamo models for planets other than Earth. Space Sci. Rev. 152, 617–649.
  * Stellmach et al. [2014] Stellmach, S., Lischper, M., Julien, K., Vasil, G., Cheng, J. S., Ribeiro, A., King, E. M. & Aurnou, J. M. 2014 Approaching the asymptotic regime of rapidly rotating convection: boundary layers versus interior dynamics. Phys. Rev. Lett. 113, 254501.
  * Stevens et al. [2018] Stevens, R. J. A. M., Blass, A., Zhu, X., Verzicco, R. & Lohse, D. 2018 Turbulent thermal superstructures in Rayleigh-Bénard convection. Phys. Rev. Fluids  3 (4), 041501.
  * Vogt et al. [2021] Vogt, T, Horn, S & Aurnou, J M 2021 Oscillatory thermal–inertial flows in liquid metal rotating convection. J. Fluid Mech.  911, A5.
  * Yan & Calkins [2022a] Yan, M. & Calkins, M. A. 2022a Asymptotic behaviour of rotating convection-driven dynamos in the plane layer geometry. J. Fluid Mech.  951, A24.
  * Yan & Calkins [2022b] Yan, M. & Calkins, M. A. 2022b Strong large scale magnetic fields in rotating convection-driven dynamos: The important role of magnetic diffusion. Phys. Rev. Res.  4 (1), L012026.
  * Yan et al. [2021] Yan, M., Tobias, S. M. & Calkins, M. A. 2021 Scaling behaviour of small-scale dynamos driven by Rayleigh–Bénard convection. J. Fluid Mech.  915, A15.