# Parking on supercritical geometric Bienaymé–Galton–Watson trees

Linxiao Chen & Alice Contat Université Sorbonne Paris-Nord.
<EMAIL_ADDRESS>Sorbonne Paris-Nord & FSMP.
<EMAIL_ADDRESS>
###### Abstract

Consider a supercritical Bienaymé–Galton–Watson tree $\mathcal{T}$ with
geometric offspring distribution. Each vertex of this tree represents a
parking spot which can accommodate at most one car. On the top of this tree,
we add $(A_{u}:u\in\mathcal{T})$ i.i.d. non negative integers sampled
according to a given law $\mu$, which are the car arrivals on $\mathcal{T}$.
Each car tries to park on its arriving vertex and if the spot is already
occupied, it drives towards the root and takes the first available spot. If no
spot is found, then it exits the tree without parking. In this paper, we
provide a criterion to determine the phase of the parking process
(subcritical, critical, or supercritical) depending on the generating function
of $\mu$.

##  1 Introduction

The parking process was introduced by Konheim & Weiss [10] in the case of an
oriented line. Recently, its generalization on _rooted trees_ was initiated by
Lackner & Panholzer [11]. Since then, the phase transition of this model
received much attention and was proved and located on _critical_
Bienaymé–Galton–Watson trees in an increasing level of generality [9, 8, 6,
7].

In the case of _supercritical_ Bienaymé–Galton–Watson trees, the only model
for which the threshold of the transition is explicitly known is a
deterministic tree, namely that of the infinite binary tree [1]. The existence
of the phase transition is however proved in a more general context [9, 2] and
the model is closely related to Derrida–Retaux model [5]. In this paper, we
locate and study the phase transition of the parking process on supercritical
Bienaymé–Galton–Watson trees with a geometric offspring distribution using
essentially enumeration results provided by Chen [4].

#### Parking on trees.

We focus here on the case where the underlying tree is a supercritical
Bienaymé–Galton–Watson tree $\mathcal{T}$ with a geometric offspring
distribution i.e. a tree which starts from one individual and where each
individual reproduces independently according to the offspring distribution

$\nu:=\sum_{k\geqslant 0}q^{k}(1-q)\delta_{k}$

with parameter $1/2<q<1$, so that $\nu$ has mean $q/(1-q)>1$. We denote by
$|\mathcal{T}|$ its size (number of vertices). It is well-known that this tree
has a positive probability to be infinite since $q>1/2$, which is not the case
when $q\leqslant 1/2$. On the top of this tree, we add
$(A_{u}:u\in\mathcal{T})$ i.i.d. non negative integers with law $\mu$, which
represent the car arrivals on $\mathcal{T}$. Each car tries to park at its
arriving vertex, but if it is already occupied, it drives towards the root
until it encounters an empty vertex and parks there. If it arrives at the root
without parking, then it exits the root and contributes to the _flux_ of
outgoing cars. To avoid cases where all cars park on their arrival node, we
assume $\mu(\\{0,1\\})<1$. When the tree is finite, the final configuration
does not depend on the order in which cars are parked. When the tree is
infinite, we circumvent the possible issue by parking the cars layers by
layers, see for example [1, Section 2.1] or the proof of Lemma 2. We introduce

$X:=\mbox{ the number of cars visiting the root of }\mathcal{T},$

so that the flux of outgoing cars equals $(X-1)_{+}=\max(X-1,0)$. In [9,
Theorem 3.4], Goldschmidt and Przykucki proved that there are two possible
regimes for the parking process on the supercritical Bienaymé–Galton–Watson
trees $\mathcal{T}$ depending on the two laws $\mu$ and $\nu$:

  * •

Either $\mathbb{E}[X]<\infty$ _(subcritical regime)_ ,

  * •

Or $X=\infty$ as soon as $\mathcal{T}$ is infinite _(supercritical regime)_.

The goal of this paper is to provide a criterion to determine if the parking
process is subcritical or supercritical depending on the parameter $q$ and the
car arrivals distribution $\mu$. To state our theorem, we need to introduce
$G$ the generating function of the car arrival law $\mu$

$G(t):=\sum_{k\geqslant 0}\mu_{k}x^{k}.$

######  Theorem 1.

Suppose that there exists $t_{c}$ such that

$t_{c}:=\inf\\{t>0,(G(t)-tG^{\prime}(t))^{2}=2t^{2}G(t)G^{\prime\prime}(t)\\}.$

Then the parking process is subcritical if and only if

$t_{c}>1\quad\mbox{ and
}\quad\frac{t_{c}G(t_{c})}{\varphi(t_{c})^{2}}\leqslant q(1-q),$

where $\varphi(y)=(y+1)G(y)-y(y-1)G^{\prime}(y)$.

Let us make a couple of remarks on this theorem. First to determine if
$t_{c}>1$, it suffices to check if $(G(1)-1G^{\prime}(1))^{2}>2\cdot
1^{2}G(1)G^{\prime\prime}(1)$, which is equivalent to $2\sigma^{2}+m^{2}<1$
where $m$ and $\sigma^{2}$ are respectively the mean and the variance of
$\mu$. Moreover, the second inequality is always satisfied for $q=1/2$. In
particular, when $q\searrow 1/2$, we recover the same criterion for the phase
transition as that proved by Curien and Hénard in [8, Theorem 1]. However, the
nature of the phase transition is different when $q=1/2$, since the tree
$\mathcal{T}$ is almost surely finite. In this case, the flux of outgoing cars
has an infinite expectation in the supercritical and critical phase of the
parking process, whereas this expectation is finite in the subcritical phase,
see also [2].

Moreover, when $t_{c}\leqslant 1$, the parking process is always
supercritical. As an example, if $\mu$ is a car arrival distribution such that
the radius of convergence of $G$ equals $1$, then for all $\alpha\in(0,1]$,
the parking process with car arrival distribution
$\mu_{\alpha}=(1-\alpha)\delta_{0}+\alpha\mu$ is supercritical. In other
words, the phase transition is trivial for the stochastically increasing
family of car arrival distributions $(\mu_{\alpha}:\alpha\in[0,1])$ since it
occurs at $\alpha_{c}:=0$: if $\alpha=0$, the parking process is subcritical
(no cars arriving) and if $\alpha\geqslant 0$, then the parking process is
supercritical.

When $t_{c}$ does not exist, we are in the _dense phase_ in the classification
of Chen [4]. Our criterion remains valid by verifying the two inequalities for
$\widetilde{t}_{c}$ the radius of convergence of $G$.

The proof follows essentially the same ideas as in [1] and rely on a
combinatorial decomposition of the final configuration of parking into
clusters of parked cars and on the enumeration of plane _fully parked trees_
by Chen [4]. We use the fact that the flux of outgoing cars is almost surely
finite in the subcritical regime to obtain an equation involving the
generating function of fully parked trees which characterizes the subcritical
regime, see Proposition 2.

#### Examples.

Let us give some examples of explicit laws where the criterion to determine
the threshold of the parking process is explicit. To make this criterion more
intelligible, it is more convenient to fix the car arrival law and determine
the range of parameter $q$ for which the parking process is subcritical. We
will come back to this examples in Section 3.2.
Binary arrivals. The simplest example of non-trivial car arrival distribution
is the binary distribution $\mu_{\alpha}$: each vertex receives independently
either two cars with probability $\alpha/2$ or no car with probability
$1-\alpha/2$ for some parameter $\alpha\in[0,1]$. In this case, the generating
function of $\mu_{\alpha}$ is $G(t):=(1-\alpha/2)+\alpha/2t^{2}$. Then, the
parking process is subcritical if and only if

$\alpha\leqslant 2-\sqrt{3}\quad\mbox{ and }\quad
q\leqslant\frac{1}{2}\left(1+\sqrt{1-\frac{6\sqrt{2\sqrt{3}+3}}{\sqrt{\alpha(2-\alpha)}\left(3+\sqrt{2\sqrt{3}-3}\sqrt{\frac{2-\alpha}{\alpha}}\right)^{2}}}\right).$

Geometric arrivals. If $\mu_{\alpha}=\sum_{k\geqslant 0}p^{k}(1-p)\delta_{k}$
is a geometric distribution with mean $\alpha=p/(1-p)$, then the parking
process is subcritical if and only if

$\alpha\leqslant\frac{1}{3}\quad\mbox{ and }\quad
q\leqslant\frac{1}{2}\left(1+\frac{(1-3\alpha)^{3/2}}{1+9\alpha}\right).$

Poisson arrivals. Lastly, if $\mu_{\alpha}$ is Poisson with mean $\alpha$,
then the parking process is subcritical if and only if

$\alpha\leqslant\sqrt{2}-1\quad\mbox{ and
}q\leqslant\frac{1}{2}\left(1+\sqrt{1-\frac{2(\sqrt{2}-1)\alpha\mathrm{e}^{\alpha-(\sqrt{2}-1)}}{(\alpha+3-2\sqrt{2})^{2}}}\right).$

Non generic dilute distributions. Another class of examples is when the car
arrivals distribution is non generic. More precisely, let us take the
generating function $G$ of the form

$G(t)=P(t)+C\left(1-\frac{t}{\rho}\right)^{\alpha}.$

where $P$ is a polynomial, and we fix a negative constant $C<0$ , a radius of
convergence $\rho>1$ and an exponent $\alpha\in(2,3)$. We assume that $G$ is
the generating function of a critical law, in the sense that $t_{c}=\rho$ and
the parameter $q$ of the tree $\mathcal{T}$ is such that the inequality on the
right in Theorem 1 is an equality. Then, the generating function of $X$ the
number of cars visiting the root has a radius of convergence $t_{c}$. In
particular, the tail of the law of $X$ decays exponentially fast. However,
this is not the case for the law of the size $|\mathcal{C}(\varnothing)|$ of
the cluster of the root. For example, using the asymptotic provided by Chen in
[4, Theorem 1], we can show that conditionally on $X=1$, the probability that
$|\mathcal{C}(\varnothing|)$ has size $n$ decays polynomially:

$\mathbb{P}\left(|\mathcal{C}(\varnothing)|=n|X=1\right)\underset{n\to\infty}{\sim}\mathrm{cst}\cdot
n^{\frac{-2\alpha+1}{\alpha-1}},$

for some constant $\mathrm{cst}>0$ that depends on $\alpha$, see Section 3.2.
Note that the exponent $(-2\alpha+1)/(\alpha-1)$ is increasing from $-3$ to
$-5/2$ for $\alpha\in(2,3)$.

#### Outline of the paper.

In Section 2, we show the dichotomy between the two regimes for the parking
process, as well as our main decomposition which allows us to characterizes
these two regimes using the generating function of _fully parked trees_.
Section 3 is devoted to the proof of our main Theorem 1. We also detail in
this section some examples to which our theorem applies and the critical
exponents which appear in this examples.

#### Acknowledgments.

We thank Nicolas Curien for stimulating discussions. The second author
acknowledges partial support from SuPerGRandMa, the ERC CoG 101087572.

##  2 Phase transition

In this section, we prove some “probabilistic” properties of the parking
process that characterize the two different regimes of this model.

###  2.1 Two regimes

We start by proving the dichotomy between the two different regimes of the
parking process. This distinction already appears in [9, Theorem 3.4] but in
our case, we show here a slightly more narrow property to characterize the
subcritical regime.

######  Lemma 1 (Dichotomy subcritical/supercritical regime).

We have the following dichotomy :
Subcritical regime. Either
$\left(\frac{q}{1-q}\right)^{n}\mathbb{P}\left(X\geqslant n\right)$ is
bounded,
Supercritical regime. Or $X=\infty$ as long as $\mathcal{T}$ is infinite.

Notice in particular that the number of cars visiting the root is almost
surely finite and has a finite expectation in the subcritical regime.

###### Proof.

Let us denote by $Z_{n}$ the number of vertices at level $n$ (at distance $n$
from the root) in the tree $\mathcal{T}$. Since the offspring distribution
$\nu$ of $\mathcal{T}$ has mean $q/(1-q)$ and has bounded variance, then
$(q/(1-q))^{-n}Z_{n}$ is a positive martingale which is bounded in $L^{2}$.
Thus, it converges is $L^{2}$ and the probability that the limit equals $0$ is
the same as the probability that the tree $\mathcal{T}$ is finite. As a
consequence, we have

$\lim_{\varepsilon\to
0}\lim_{n\to\infty}\mathbb{P}\left(Z_{n}\geqslant\varepsilon\left(\frac{q}{1-q}\right)^{n}\right)=\mathbb{P}\left(|\mathcal{T}|=\infty\right).$

Assume now that the sequence $(q/(1-q))^{n}\mathbb{P}\left(X\geqslant
n\right)$ is not bounded. Then this is also the case for the sequence
$(q/(1-q))^{n}\mathbb{P}\left(X\geqslant n+k\right)$ for any integer $k$. Then
on each of the $Z_{n}$ vertices at level $n$, independent copies of
Bienaymé–Galton–Watson trees are attached and for each vertex $u$ at height
$n$, we can consider $X(u)$ the random number of cars visiting the vertex $u$.
We obtain the upper bound

$\displaystyle\mathbb{P}\left(\bigcap_{u\text{ at height
}n}\left\\{X(u)\leqslant n+k\right\\}\mbox{ and
}Z_{n}\geqslant\varepsilon\left(\frac{q}{1-q}\right)^{n}\right)$
$\displaystyle\leqslant$
$\displaystyle(1-\mathbb{P}\left(X>n+k\right))^{\varepsilon\left(\frac{q}{1-q}\right)^{n}}$
$\displaystyle\leqslant$
$\displaystyle\exp\left(-\varepsilon\left(\frac{q}{1-q}\right)^{n}\mathbb{P}\left(X>n+k\right)\right).$

and the right-hand side goes to $0$ along a subsequence. But on the complement
of the event on the left-hand side, either $X\geqslant k$ or
$Z_{n}\leqslant\varepsilon({q}/{(1-q)})^{n}$. Thus, with probability at least
$\mathbb{P}\left(|\mathcal{T}|=\infty\right)$, we have $X\geqslant k$. Since
$k$ is arbitrary, this implies that $X=\infty$ with probability at least
$\mathbb{P}\left(|\mathcal{T}|=\infty\right)$. Moreover, the number of cars
visiting the root can not be infinite when the tree is infinite, thus
$X=\infty$ with probability at most
$\mathbb{P}\left(|\mathcal{T}|=\infty\right)$, which concludes the proof.

∎

Note that $X$ the number of cars visiting the root is at least the number of
cars arriving at the root. In particular, if the radius of convergence of $G$
the generating function of the car arrivals is strictly smaller than
$q/(1-q)$, then it is also the case for the generating function of the
variable $X$, which implies that the parking process is supercritical.

We are now interested in the possible sizes of the clusters that we can
observe depending on the regime. We call _black_ the clusters of parked
vertices and _white_ the clusters of empty vertices. We now prove that all
black clusters are finite in the subcritical regime (which is not true for all
white clusters). This is the analog of [1, Proposition 1] in our case and the
proof essentially follows the same line. We adapt it in our case to keep our
paper self-contained.

######  Lemma 2.

In the subcritical regime, there is no infinite black cluster.

###### Proof.

This proof is mostly an adaptation of the proof of [1, Proposition 1]. As
mentioned in [1, Section 2.1], we park the lowest cars first to avoid
Abelianity issues when the tree $\mathcal{T}$ is infinite. More precisely, for
each $n\geqslant 0$, we denote by $[\mathcal{T}]_{n}$ the finite tree obtained
by cutting $\mathcal{T}$ above height $n$, so that $[\mathcal{T}]_{n}$ is made
of all vertices of $\mathcal{T}$ of height smaller or equal than $n$. Then,
parking the cars arriving on $[\mathcal{T}]_{n}$, we can define the variables

$X_{n}(u)\quad u\in[\mathcal{T}]_{n}$

denoting the number of cars visiting the vertex $u\in[\mathcal{T}]_{n}$ in the
final configuration of parking in $[\mathcal{T}]_{n}$. Since for all
$u\in[\mathcal{T}]_{n}$, the sequence $(X_{n}(u):n\geqslant 0)$ is non-
decreasing, we can define

$X(u):=\lim_{n\to\infty}X_{n}(u)$

which represents the number of cars visiting $u\in\mathcal{T}$ even when
$\mathcal{T}$ is infinite.

Now the parking rules are clear even on infinite trees, we assume the tree
$\mathcal{T}$ is infinite (otherwise they can not be infinite clusters) and
the law $\mu$ is subcritical. Thus, all variables $(X(u):u\in\mathcal{T})$ are
almost surely finite and we denote by $\mathcal{C}(\varnothing)$ the cluster
of the root of $\mathcal{T}$. Let us fix $p\geqslant 0$ and consider
$\mathcal{E}=\\{X(\varnothing)=p\mbox{ and }\mathcal{C}(\varnothing)\mbox{ is
infinite}\\}$. We can construct a sequence of stopping times
$\theta_{1}<\theta_{1}<\dots$ such that $\theta_{1}:=\inf\\{n\geqslant
0:X_{n}(\varnothing)=p\\}$ and for all $i\geqslant 1$, we define sequentially
$\theta_{i+1}:=\infty\\{n>\theta_{i}:\varnothing\leftrightarrow\partial[\mathcal{T}]_{n}\\}$
meaning the $\theta_{i+1}$ is the first $n>\theta_{i}$ such that there exists
a path connecting the root $\varnothing$ and the level $n$ in
$[\mathcal{T}]_{n}$ which satisfies $X_{n}(u)\geqslant 1$ for all vertices $u$
in this path. Let us denote by $v_{i}$ the leftmost vertex at level
$\theta_{i}$ which is connected to the root by such a path (if there is one)
for $i\geqslant 2$. Note that on the event $\mathcal{E}$, then all stopping
times $\theta_{i}$ are finite and thus, for all $i\geqslant 2$, the vertex
$v_{i}$ is well-defined and its children have at most $1$ car arriving on it.
In particular, we obtain for all $n\geqslant 2$

$\displaystyle\mathbb{P}\left(\mathcal{E}\right)$ $\displaystyle\leqslant$
$\displaystyle\prod_{i=2}^{n}\mathbb{P}\left(\mbox{all children of
}v_{i}\mbox{ receive at most one car}|\theta_{i}<\infty\right)$
$\displaystyle\leqslant$ $\displaystyle\left(\sum_{k\geqslant
0}q^{k}(1-q)(\mu_{0}+\mu_{1})^{k}\right)^{n-1}=\left(\frac{1-q}{1-q(\mu_{0}+\mu_{1})}\right)^{n-1}.$

Since we assumed $\mu(\\{0,1\\})<1$ and $n$ is arbitrary, the fraction this
implies $\mathbb{P}\left(\mathcal{E}\right)=0$ and this is true for all value
of $p$, which concludes the proof.

∎

###  2.2 Fully parked trees

This section presents a description of the fully parked components and mainly
contains the results from Section 2 in [4].

Recall our notation $(a_{x}:x\in\mathbf{t})$ for the car arrivals on the
vertices of a tree $\mathbf{t}$. A _fully parked tree_ is a tree decorated
with its car arrivals configuration such that in the final configuration, all
vertices are occupied spots and they may be outgoing cars. We then denote by
$\mathbb{T}_{n}^{(k)}$ the set of (plane) fully parked trees with $n$ vertices
and $k$ outgoing cars (hence $n+k$ arriving cars) and introduce the bivariate
generating function

$F(x,y)=\sum_{n\geqslant 1}\sum_{k\geqslant
0}\sum_{\mathbf{t}\in\mathbb{T}_{n}^{(k)}}w(\mathbf{t})x^{n}y^{k}\quad\mbox{
where}\quad w(\mathbf{t}):=\prod_{x\in\mathbf{t}}\mu_{a_{x}}.$

The first important tool is a recursive equation obtained from an “à la Tutte”
recursive decomposition, i.e. by decomposing each fully parked tree according
to the splitting at its root. Indeed, if the root of
$\mathbf{t}\in\mathbb{T}_{n}^{(k)}$ has $j\geqslant 0$ children, then all
these $j$ children are the respective roots of a sequence
$(\mathbf{t}_{i}:1\leqslant i\leqslant j)$ of fully parked trees. Moreover, if
for all $1\leqslant i\leqslant j$, the tree $\mathbf{t}_{i}$ belongs to
$\mathbb{T}_{n_{i}}^{(k_{i})}$ then the fact that $\mathbf{t}$ has $n$
vertices and $k$ outgoing cars implies that $n=1+\sum_{i=1}^{j}n_{i}$ and
$k=\sum_{i=1}^{j}k_{i}+a_{\varnothing}-1$ where $a_{\varnothing}$ is the
number of cars arriving at the root of $\mathbf{t}$, see Figure 1. Note that
since the root of $\mathbf{t}$ is an occupied spot, then we can not have
$a_{\varnothing}=0$ and $k_{i}=0$ for all $1\leqslant i\leqslant j$
simultaneously. Translating this observation on $F$ gives the following
equation (see also [4, Equation 21])

$F(x,y)=\frac{x}{y}\left(\frac{G(y)}{1-F(x,y)}-\frac{G(0)}{1-F(x,0)}\right).$
(1)

Figure 1: Illustration of the “a la Tutte” recursive decomposition. On the
left, a fully parked tree of size $n$ and $p$ outgoing cars. On the right, the
different possibility depending on the number $j$ of the children of the root
vertex.

Note that this equation characterizes $F$ since the coefficient in front of
$x^{n}$ on the left-hand side only depends on the coefficients of $F(x,y)$ in
front of $x^{k}$ for $0\leqslant k\leqslant n-1$. It involves not only
$F(x,y)$ but also $F(x,0)$ its specialization at $y=0$. We thus apply the
usual _kernel method_ [3] to solve this equation as it is done in [4, Section
2]. We introduce
$P(f,f_{0},x,y)=\frac{G(y)}{1-f}-\frac{G(0)}{1-f_{0}}+\frac{xf}{y}$, so that
Equation (1) is equivalent to $P(F(x,y),F(x,0),x,y)=0$. The key idea is to
find a function $Y=Y(x)$ such that $\partial_{f}P(F(x,Y(x)),F(x,0),x,Y(x))=0$
so that differentiating $P(F(x,Y(x)),F(x,0),x,Y(x))=0$ with respect to $y$, we
also obtain that $\partial_{y}P(F(x,Y(x)),F(x,0),x,Y(x))=0$. Using the
explicit form for the function $P$, we then have

$\frac{G(Y)}{1-F(x,Y)}=\frac{Y}{x},\qquad\frac{G^{\prime}(Y)}{1-F(x,Y)}=\frac{F(x,Y)}{x},\qquad\mbox{and}\qquad\frac{G(Y)}{1-F(x,Y)}-\frac{G(0)}{1-F(x,0)}=\frac{YF(x,Y)}{x}.$
(2)

The existence of a solution for $Y$ is guaranteed by [3] and eliminating the
variable $F(x,Y)$, we can express $x$ and $F(x,0)$ in terms of $Y$ and get

$x=\frac{YG(Y)}{(G(Y)+YG^{\prime}(Y))^{2}}\qquad\mbox{and}\qquad
F(x,0)=1-\frac{G(0)G(Y)}{(G(Y)-YG^{\prime}(Y))^{2}}.$ (3)

Plugging this back into Tutte’s equation (1), we can obtain a parametrization
of $F(x,y)$. In particular, we can obtain its value and radius of convergence
in $x$ at $y=1$ , which will be useful thereafter and which we sum up in the
following proposition. See also [4, Proposition 1].

######  Proposition 1.

As in Theorem 1, we assume that there exists $t_{c}$ such that

$t_{c}:=\inf\\{t>0,(G(t)-tG^{\prime}(t))^{2}=2t^{2}G(t)G^{\prime\prime}(t)\\}.$

We also suppose that $t_{c}>1$. Then, writing
$x=\hat{x}(Y)=\frac{YG(Y)}{(G(Y)+YG^{\prime}(Y))^{2}},$ we have

$\displaystyle
F(x,1)=\hat{F}(Y,1)=\left\\{\begin{array}[]{lll}1+\frac{\varphi(Y)}{2(G(Y)+YG^{\prime}(Y))}\left(-\sqrt{1-\frac{4YG(Y)}{\varphi(Y)^{2}}}-1\right)&\text{
when }&0\leqslant Y\leqslant 1,\\\
1+\frac{\varphi(Y)}{2(G(Y)+YG^{\prime}(Y))}\left(\sqrt{1-\frac{4YG(Y)}{\varphi(Y)^{2}}}-1\right)&\text{
when }&1<Y\leqslant t_{c},\end{array}\right.$

where $\varphi(y)=(y+1)G(y)-y(y-1)G^{\prime}(y)$. Moreover, the radius of
convergence of $x\mapsto F(x,1)$ is
$\hat{x}(t_{c})=\frac{t_{c}G(t_{c})}{(G(t_{c})+t_{c}G^{\prime}(t_{c}))^{2}}$.

###### Proof.

The parametrization follows easily from the kernel method applied to Tutte’s
equation (1). In particular, it is well-defined and the quantity inside the
square root is non negative for all $0\leqslant Y\leqslant t_{c}$. To see that
$\hat{x}(t_{c})$ is the radius of convergence of $F(x,1)$, we first observe
that the parametrization implies that this radius is at least
$\hat{x}(t_{c})$, see [3]. If $G$ or one of its first two derivative of $G$ is
infinite at $t_{c}$, then the result is clearly true. Otherwise, thanks to our
parametrization, we can show that the function $F(x,1)$ and its first
derivative have a finite limit as $x\to\hat{x}(t_{c})$, see also [4, Lemma
14]. Moreover, by tedious computations, we can show that

$\lim_{x\to\hat{x}(t_{c})}\frac{\partial^{2}F(x,1)}{\partial x^{2}}=\lim_{y\to
t_{c}}\frac{(G(y)+yG^{\prime}(y))^{5}}{yG(y)}\left(\sqrt{1+\frac{4yG(y)}{\varphi(y)^{2}-4yG(y)}}-1\right)\cdot\frac{1}{(G(y)-yG^{\prime}(y))^{2}-2y^{2}G(y)G^{\prime\prime}(y)}$

The first factor has a finite limit whereas the denominator vanishes on the
right-hand side. This concludes the proof.

∎

###  2.3 Decomposition into fully-parked components

In this section, we present the main decomposition from which we derive our
characterization of the subcritical regime for the parking process. We
remember that $X$ is the number of cars visiting the root of $\mathcal{T}$
with $p_{\circ}=\mathbb{P}\left(X=0\right)$ and we write
$\mathcal{C}(\varnothing)$ for the cluster of the root (either white (void) or
black (parked)) in $\mathcal{T}$ after parking.

Let us fix $k\geqslant 0$. For all $x\geqslant 0$, we denote by
$[y^{k}]F(x,y)$ denote the coefficient in front of $y^{k}$ in the (formal)
power series $y\mapsto F(x,y)$. For all $x\geqslant 0$ such that
$[y^{k}]F(x,y)<\infty$, we say that $\mathcal{F}$ is a (random) fully parked
tree with $k$ outgoing cars and $x$-Boltzmann law if for all $n\geqslant 0$
and $\mathbf{t}\in\mathbb{T}_{n}^{(k)}$, we have

$\mathbb{P}\left(\mathcal{F}=\mathbf{t}\right)=\frac{1}{[y^{k}]F(x,y)}x^{n}w(\mathbf{t}).$

In the subcritical regime, we can link the variable $X$ and the coefficients
of $F$ through the following lemma.

######  Lemma 3.

When the law $\mu$ is subcritical for the parking process, for all $k\geqslant
0$, we have

$\mathbb{P}\left(X=k+1\right)=\frac{1-qp_{\circ}}{q}[y^{k}]F\left(\frac{q(1-q)}{(1-qp_{\circ})^{2}},y\right).$

Moreover, conditionally on $X=k+1$, the cluster of the root
$\mathcal{C}(\varnothing)$ is a fully parked tree with $k$ outgoing cars and
$\frac{q(1-q)}{(1-qp_{\circ})^{2}}$-Bolztmann law.

###### Proof.

The proof follows from a simple decomposition according to the cluster of the
root. We emphasize here the important role of the geometric offspring
distribution here: in the tree $\mathcal{T}$, each vertex has a probability
$q^{k}(1-q)$ to have $k$ children, which is equivalent to saying that each
edge carries a weight $q$ and each vertex carries a weight $(1-q)$. Let us fix
a finite fully parked tree $\mathbf{t}\in\mathbb{T}_{n}^{(k)}$ and compute the
probability that $\mathcal{C}(\varnothing)=\mathbf{t}$. The weight of this
tree is the product of the weight $w(\mathbf{t})$ of the car arrivals on it,
the weight $(1-q)^{n}$ of its $n$ vertices and $q^{n-1}$ for its $n-1$ edges.
Moreover, since $\mathbf{t}$ has $n$ vertices, there are $2n-1$ corners where
there can be empty spots. Each empty spot comes with a weight $q$ for edge
between the vertex and a vertex of $\mathbf{t}$ multiplied by $p_{\circ}$
since it can be seen as the root of a Bienaymé–Galton–Watson with an empty
root. The weight of each corner is then $\frac{1}{1-qp_{\circ}}$. See Figure 2
for an example. We thus obtain

$\mathbb{P}\left(\mathcal{C}(\varnothing)=\mathbf{t}\right)=w(\mathbf{t})(1-q)^{n}q^{n-1}\left(\frac{1}{1-qp_{\circ}}\right)^{2n-1}=\frac{(1-qp_{\circ})}{q}\cdot\left(\frac{q(1-q)}{(1-qp_{\circ})^{2}}\right)^{n}w(\mathbf{t})$

Figure 2:  On the left, a tree together with its car arrivals configuration.
On the right, the (black) cluster of root of this tree together with all empty
spots attached to it in the initial tree. All this empty spots are the root of
independent Bienaymé–Galton–Watson trees where in the final configuration, the
root is empty. This cluster is a fully parked tree with $7$ vertices and there
are $8$ empty spots attached to it. The probability to observe this cluster is
$(\mu_{0}^{2}\mu_{1}^{2}\mu_{2}^{2}\mu_{3})q^{6}(1-q)^{7}(qp_{\circ})^{8}$

Moreover, we know from Lemma 2 that all black clusters are finite in the
subcritical regime. Thus, for all $k\geqslant 0$

$\displaystyle\mathbb{P}\left(X=k+1\right)$ $\displaystyle=$
$\displaystyle\mathbb{P}\left(X=k+1\mbox{ and
}|\mathcal{C}(\varnothing)|<\infty\right)$ $\displaystyle=$
$\displaystyle\sum_{n\geqslant
1}\sum_{\mathbf{t}\in\mathbb{T}_{n}^{(k)}}\mathbb{P}\left(|\mathcal{C}(\varnothing)|=\mathbf{t}\right)$
$\displaystyle=$ $\displaystyle\sum_{n\geqslant
1}\sum_{\mathbf{t}\in\mathbb{T}_{n}^{(k)}}\frac{(1-qp_{\circ})}{q}\cdot\left(\frac{q(1-q)}{(1-qp_{\circ})^{2}}\right)^{n}w(\mathbf{t})=\frac{(1-qp_{\circ})}{q}[y^{k}]F\left(\frac{q(1-q)}{(1-qp_{\circ})^{2}},y\right),$

which concludes the proof. ∎

Moreover, in the subcritical regime, the variable $X$ is almost surely finite,
thus splitting according to the values of $X$, we obtain

$1=p_{\circ}+\sum_{k\geqslant
0}\mathbb{P}(X=k+1)=p_{\circ}+\frac{1-qp_{\circ}}{q}F\left(\frac{q(1-q)}{(1-qp_{\circ})^{2}},1\right).$
(5)

Note that the function $p\mapsto
p+\frac{1-qp}{q}F\left(\frac{q(1-q)}{(1-qp)^{2}},1\right)$ is an increasing
function of $p$. Thus the solution to this equation is unique if there is one.
This gives us the following characterization of the subcritical regime.

######  Proposition 2 (F-characterization of the subcritical regime).

The law $\mu$ is subcritical for the parking process on a
Bienaymé–Galton–Watson tree with geometric offspring distribution with
parameter $q$ if and only if there exists a positive solution $p_{\circ}>0$ to
the equation

$\frac{1-qp}{q}\cdot F\left(\frac{q(1-q)}{(1-qp)^{2}},1\right)+p=1.$ (6)

The proof follows the same line as [1, Proposition 2]

###### Proof.

If the parking process is subcritical, then $p_{\circ}$ is solution of
Equation (6) by the above computations. Reciprocally suppose that there exists
a solution $p_{\circ}$ to the equation (6). Then there exist a random variable
$Z$ whose generating function is

$\mathbf{Z}(y):=\frac{1-qp_{\circ}}{q}\cdot
F\left(\frac{q(1-q)}{(1-qp_{\circ})^{2}},y\right)+p_{\circ}.$

Then, writing $x_{\circ}:=q(1-q)/(1-qp_{\circ})^{2}$ and $F_{0}(x):=F(x,0)$,
we compute

$\displaystyle\frac{1}{y}\left(\frac{(1-q)G(y)}{\left(1-q\mathbf{Z}(y)\right)}-\frac{(1-q)G(0)}{\left(1-q\mathbf{Z}(0)\right)}\right)+\frac{(1-q)G(0)}{\left(1-q\mathbf{Z}(0)\right)}$
$\displaystyle=$
$\displaystyle\frac{1}{y}\left(\frac{(1-q)G(y)}{(1-qp_{\circ})\left(1-F\left(x_{\circ},y\right)\right)}-\frac{(1-q)G(0)}{(1-qp_{\circ})\left(1-F\left(x_{\circ},0\right)\right)}\right)+\frac{(1-q)G(0)}{(1-qp_{\circ})\left(1-F\left(x_{\circ},0\right)\right)}$
$\displaystyle\overset{\eqref{eq:tutte}}{=}$
$\displaystyle\frac{1-q}{1-qp_{\circ}}\cdot\frac{1}{x_{\circ}}F(x_{\circ},y)+\frac{(1-q)G(0)}{(1-qp_{\circ})\left(1-F\left(x_{\circ},0\right)\right)}$
(7)

Moreover, from Equation (6), we obtain

$F\left(x_{\circ},1\right)=\frac{(1-p_{\circ})q}{1-p_{\circ}q}\quad\mbox{and}\quad
1-F\left(x_{\circ},1\right)=\frac{1-q}{1-p_{\circ}q}.$

Thus, using Tutte’s equation (1) with $y=1$, we obtain

$\displaystyle\frac{(1-q)G(0)}{(1-qp_{\circ})\left(1-F\left(x_{\circ},0\right)\right)}$
$\displaystyle=$
$\displaystyle\frac{1-q}{1-qp_{\circ}}\left(\frac{1}{1-F\left(x_{\circ},1\right)}-\frac{F(x_{\circ},1)}{x_{\circ}}\right)$
$\displaystyle=$ $\displaystyle
1-\frac{(1-q)}{(1-qp_{\circ})}\cdot\frac{(1-qp_{\circ})^{2}}{q(1-q)}\cdot\frac{(1-p_{\circ})q}{(1-p_{\circ}q)}$
$\displaystyle=$ $\displaystyle 1-(1-p_{\circ})=p_{\circ}.$

Hence, plugging this equality in (7), we obtain

$\mathbf{Z}(y)=\frac{1}{y}\left(\frac{(1-q)G(y)}{\left(1-q\mathbf{Z}(y)\right)}-\frac{(1-q)G(0)}{\left(1-q\mathbf{Z}(0)\right)}\right)+\frac{(1-q)G(0)}{\left(1-q\mathbf{Z}(0)\right)}.$

This identity is equivalent to the following recursive distributional equation
for $Z$:

$Z\overset{(d)}{=}\left(\sum_{i=1}^{Y}Z_{i}+A-1\right)_{+},$

where on the right-hand side, the variables $(Z_{i}:i\geqslant 1)$ are i.i.d.
copies of $Z$, the variable $Y$ has law $\nu$ the offspring distribution of
the tree, the variable $A$ has law $\mu$ the car arrivals distribution and all
variables are independent.

This recursive equation enables us to construct the tree $\mathcal{T}$ (using
the variable $Y$) together with its car arrivals (using the variable $A$)
layers by layers in a coherent matter ans so that in the end, the flux of out
cars has law $Z$. In particular, it implies that the parking process with car
arrival law $\mu$ on the tree $\mathcal{T}$ is subcritical. ∎

##  3 Probabilistic consequences

###  3.1 Location of the thereshold: Theorem 1

We have now all the tools to prove Theorem 1. Recall that thanks to
Proposition 2, the parking process is subcritical if and only if (6) has a
positive solution $p_{\circ}$. Let us write
$x_{\circ}=q(1-q)/(1-qp_{\circ})^{2}$, Equation (6) is equivalent to

$\sqrt{\frac{1-q}{qx_{\circ}}}(F(x_{\circ},1)-1)=-\frac{1-q}{q}.$ (8)

Note that the function $x\mapsto\sqrt{(1-q)/(qx)}(F(x,1)-1)$ is increasing in
$x$ and assuming $t_{c}\geqslant 1$, let us compute its value at $\hat{x}(1)$
using the parametrization given in Proposition 1. Recalling
$\varphi(y)=(y+1)G(y)-y(y-1)G^{\prime}(y)$, we can compute

$\displaystyle\sqrt{(1-q)/(q\hat{x}(1))}(F(\hat{x}(1),1)-1)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{(1-q)(G(1)+1G^{\prime}(1))^{2}}{qG(1)}}\cdot\frac{-\varphi(1)}{2(G(1)+1G^{\prime}(1))}$
$\displaystyle=$ $\displaystyle-\sqrt{\frac{1-q}{q}}\leqslant-\frac{1-q}{q}$

for $q\in(1/2,1)$. Thus, is $x_{\circ}$ is a solution of Equation (8), then
$x_{\circ}$ can be written as $\hat{x}(Y_{\circ})$ for some
$Y_{\circ}\geqslant 1$. In particular, the parking process is always
supercritical if $t_{c}<1$. Suppose now $t_{c}\geqslant 1$. It suffices to
check if $\sqrt{(1-q)/(qx_{c})}(F(x_{c},1)-1)\geqslant-(1-q)/q$ where $x_{c}$
is the radius of convergence of $x\mapsto F(x,1)$. Moreover, again by
Proposition 1, this radius is $\hat{x}(t_{c})$ and we obtain

$\sqrt{\frac{(1-q)\varphi(t_{c})^{2}}{q\cdot
4t_{c}G(t_{c})}}\left(\sqrt{1-\frac{4t_{c}G(t_{c})}{\varphi(t_{c})^{2}}}-1\right)\geqslant-\frac{1-q}{q}\quad\Leftrightarrow\quad\sqrt{1-\frac{4t_{c}G(t_{c})}{\varphi(t_{c})^{2}}}\geqslant
1-\sqrt{\frac{1-q}{q}}\sqrt{\frac{4t_{c}G(t_{c})}{\varphi(t_{c})^{2}}}.$

Since both sides of the inequality on the right are positive (this is a
consequence of the parametrization for the right-hand-side) and we can take
the square of both sides and substracting $1$, we obtain

$-\frac{4t_{c}G(t_{c})}{\varphi(t_{c})^{2}}\geqslant\frac{1-q}{q}\cdot\frac{4t_{c}G(t_{c})}{\varphi(t_{c})^{2}}-2\sqrt{\frac{1-q}{q}}\sqrt{\frac{4t_{c}G(t_{c})}{\varphi(t_{c})^{2}}}$

Finally, we obtain that the existence of a positive solution to Equation (6)
is equivalent to

$\frac{t_{c}G(t_{c})}{\varphi(t_{c})^{2}}\leqslant q(1-q),$

and Theorem 1 follows.

###  3.2 Examples.

We now come back to the three natural examples of car arrivals law mentioned
in the introduction, namely binary, Poisson and geometric car arrivals. Recall
that to check it $t_{c}\geqslant 1$, it suffices to verify if
$(1-G^{\prime}(1))^{2}\geqslant 2G^{\prime\prime}(1)$.
Binary arrivals. Let us start by the simplest example of non trivial car
arrival distribution, namely the binary distribution $\mu_{\alpha}$: each
vertex receives either two cars with probability $\alpha/2$ or no car with
probability $1-\alpha/2$ independently for some parameter $\alpha\in[0,1]$. In
this case, the generating function of $\mu_{\alpha}$ is
$G(t):=(1-\alpha/2)+\alpha/2t^{2}$.

$\displaystyle t_{c}\geqslant 1\ \Leftrightarrow\
(1-G^{\prime}(1))^{2}\geqslant 2G^{\prime\prime}(1)\ \Leftrightarrow\
(1-\alpha)^{2}\geqslant 2\alpha\ \Leftrightarrow\
\alpha\notin(2-\sqrt{3},2+\sqrt{3})$

Moreover, the equation $(G(t)-tG^{\prime}(t))^{2}=2G(t)G^{\prime\prime}(t)$ is
a quadratic equation in the variable $t^{2}$ and has only one solution of
$t_{c}=\sqrt{\frac{2}{\sqrt{3}-1}}\sqrt{\frac{2-\alpha}{\alpha}}$ which is
larger than 1 when $\alpha\leqslant 2-\sqrt{3}$. At this $t_{c}$ and recalling
that $\varphi(y)=(y+1)G(y)-y(y-1)G^{\prime}(y)$, we obtain

$\displaystyle\frac{t_{c}G(t_{c})}{\varphi(t_{c})^{2}}=\frac{3\sqrt{3+2\sqrt{3}}}{2\left(3+2\sqrt{2\sqrt{3}-3}\sqrt{\frac{2-\alpha}{\alpha}}\right)^{2}\sqrt{\alpha(2-\alpha)}}\leqslant\frac{1}{4}\quad\mbox{
when }\alpha\leqslant 2-\sqrt{3}$

Thus, using Theorem 1, we obtain that the parking process is subcritical if
and only if

$\alpha\leqslant 2-\sqrt{3}\quad\mbox{ and }\quad
q\leqslant\frac{1}{2}\left(1+\sqrt{1-\frac{6\sqrt{2\sqrt{3}+3}}{\sqrt{\alpha(2-\alpha)}\left(3+\sqrt{2\sqrt{3}-3}\sqrt{\frac{2-\alpha}{\alpha}}\right)^{2}}}\right).$

Geometric arrivals. Now assume that the car arrival
distribution$\mu_{\alpha}=\sum_{k\geqslant 0}p^{k}(1-p)\delta_{k}$ is a
geometric distribution with mean $\alpha=p/(1-p)$. Its generating function is
then $G(t):=1/(1+\alpha-\alpha t)$. Thus

$\displaystyle t_{c}\geqslant 1\ \Leftrightarrow\
(1-G^{\prime}(1))^{2}\geqslant 2G^{\prime\prime}(1)\ \Leftrightarrow\
(1-\alpha)^{2}\geqslant 4\alpha\ \Leftrightarrow\ \alpha\notin(1/3,1)$

Moreover, the quantity $(G(t)-tG^{\prime}(t))^{2}-2G(t)G^{\prime\prime}(t)$ is
a rational fraction in $t$ and its numerator is a linear function of $t$. We
can thus determine, for $\alpha\leqslant 1/3$,

$t_{c}=\frac{1+\alpha}{4\alpha},$

and we obtain

$\displaystyle\frac{t_{c}G(t_{c})}{\varphi(t_{c})^{2}}=\frac{27\alpha(1+\alpha)^{2}}{4(1+9\alpha)^{2}}.$

Using Theorem 1, the parking process is subcritical if and only if

$\alpha\leqslant\frac{1}{3}\quad\mbox{ and }\quad
q\leqslant\frac{1}{2}\left(1+\frac{(1-3\alpha)^{3/2}}{1+9\alpha}\right).$

Poisson arrivals. Lastly, we assume that the car arrivals distribution
$\mu_{\alpha}$ is Poisson with mean $\alpha$, which means that its generating
function is $G(t)=\exp\left(a(t-1)\right)$. Then

$\displaystyle t_{c}\geqslant 1\ \Leftrightarrow\
(1-G^{\prime}(1))^{2}\geqslant 2G^{\prime\prime}(1)\ \Leftrightarrow\
(1-\alpha)^{2}\geqslant 2\alpha^{2}\ \Leftrightarrow\
\alpha\in[-\sqrt{2}-1,\sqrt{2}-1]$

Moreover, the quantity $(G(t)-tG^{\prime}(t))^{2}-2G(t)G^{\prime\prime}(t)$ is
a quadratic polynomial in $t$ multiplied by $\mathrm{e}^{2a(t-1)}$. It cancels
twice, one of the two solutions is positive and corresponds to

$t_{c}=\frac{\sqrt{2}-1}{\alpha}.$

At this $t_{c}$, we obtain

$\frac{t_{c}G(t_{c})}{\varphi(t_{c})^{2}}=\frac{(\sqrt{2}-1)\alpha\mathrm{e}^{\alpha-(\sqrt{2}-1)}}{(\sqrt{2}a-4+3\sqrt{2})^{2}}.$

and thanks to Theorem 1, the parking process is subcritical if and only if

$\alpha\leqslant\sqrt{2}-1\quad\mbox{ and
}q\leqslant\frac{1}{2}\left(1+\sqrt{1-\frac{2(\sqrt{2}-1)\alpha\mathrm{e}^{\alpha-(\sqrt{2}-1)}}{(\alpha+3-2\sqrt{2})^{2}}}\right).$

Figure 3: From left to right, the critical value for $q$ as a function of
$\alpha$ for car arrivals with respectively binary, geometric and Poisson
distribution. Note that for the geometric and Poisson distribution, this
function has a finite derivative as $\alpha$ goes to $0$, whereas it is of the
form $1-\mathrm{cst}\sqrt{a}$ for some $\mathrm{cst}>0$ in the binary case

Stable cases. A last important example is when the car arrivals distribution
is non-generic. More precisely, we assume here that there exists a polynomial
$P$, three reals $C<0$, $\rho>1$ and $\alpha\in(2,3)$ such that

$G(t)=P(t)+C\left(1-\frac{t}{\rho}\right)^{\alpha},$

and such that the parameter $t_{c}$ in Theorem 1 is equal to $\rho$. Assume
that $G$ is the generating function of a critical law, in the sense that the
inequality on the right in Theorem 1 is an equality. Then, by Lemma 3, the
generating function of $X$ the number of cars visiting the root is
proportional to the function $y\mapsto F(x_{c},y)$ where $x_{c}$ is
$\hat{x}(t_{c})$. Using the asymptotics of Chen [4, Theorem 3], this function
has a radius of convergence $t_{c}$. In particular, the tail of the law of $X$
decays exponentially fast. However, this is not the case for the law of the
size $|\mathcal{C}(\varnothing)|$ of the cluster of the root. For example, by
Lemma 3 and conditionally on $X=1$, the generating function of the size of the
cluster of the root is $x\mapsto F(x_{c}x,0).$ In particular, using again the
asymptotic provided by Chen in [4, Theorem 3], we have the coefficient in
front of $x^{n}$ in $F(x_{c}x,0)$ is asymptotically of order

$[x^{n}]F(x_{c}x,0)\underset{n\to\infty}{\sim}C\cdot x_{c}^{n}\cdot
x_{c}^{-n}n^{\frac{-2\alpha+1}{\alpha-1}}=C\cdot
n^{\frac{-2\alpha+1}{\alpha-1}},$

for some constant $C>0$ that depends on $\alpha$. Note that this exponent
$(-2\alpha+1)/(\alpha-1)$ is increasing from $-3$ to $-5/2$ for
$\alpha\in(2,3)$.

###  3.3 Perspectives

In this work, we only consider the parking process on supercritical
Bienaymé–Galton–Watson trees with a geometric offspring distribution. However,
we believe that our technique can be extended to more general trees. The
challenge in more general cases will be to enumerate the fully parked trees
and include in this enumeration the free spots in the underlying tree. We will
thus at least need an additional catalytic variable counting the number of
adjacent vertices of the fully parked tree inside the underlying tree. We can
also imagine to include in the enumeration the possibility for a vertex to
carry more the one car. Moreover, we may make the car arrival law depend on
the degree of the vertex (in the underlying tree). For example, by allowing
car arrivals only on the leaves, this model is then very close to the Derrida-
Retaux model [5].

## References

  * [1] D. Aldous, A. Contat, N. Curien, and O. Hénard, Parking on the infinite binary tree, Probability Theory and Related Fields, (2023), pp. 1–24.
  * [2] R. Bahl, P. Barnet, and M. Junge, Parking on supercritical Galton-Watson trees, ALEA, 18 (2021), pp. 1801–1815.
  * [3] M. Bousquet-Mélou and A. Jehanne, Polynomial equations with one catalytic variable, algebraic series and map enumeration, J. Combin. Theory Ser. B, 96 (2006), pp. 623 – 672.
  * [4] L. Chen, Enumeration of fully parked trees, arXiv preprint arXiv:2103.15770, (2021).
  * [5] X. Chen, V. Dagard, B. Derrida, Y. Hu, M. Lifshits, and Z. Shi, The derrida–retaux conjecture on recursive models, Annals Probab., 49 (2021), pp. 637–670.
  * [6] A. Contat, Sharpness of the phase transition for parking on random trees, Random Structures & Algorithms, 61 (2022), pp. 84–100.
  * [7] , Parking on trees with a (random) given degree sequence and the frozen configuration model, arXiv preprint arXiv:2312.04472, (2023).
  * [8] N. Curien and O. Hénard, The phase transition for parking on Galton–Watson trees, Discrete Analysis, (2022).
  * [9] C. Goldschmidt and M. Przykucki, Parking on a random tree, Combinatorics, Probability and Computing, 28 (2019), pp. 23–45.
  * [10] A. G. Konheim and B. Weiss, An occupancy discipline and applications, SIAM Journal on Applied Mathematics, 14 (1966), pp. 1266–1274.
  * [11] M.-L. Lackner and A. Panholzer, Parking functions for mappings, Journal of Combinatorial Theory, Series A, 142 (2016), pp. 1 – 28.


# Acute sensitivity of global ocean circulation and heat content to eddy
energy dissipation time-scale

###### Abstract

The global ocean overturning circulation, critically dependent on the global
density stratification, plays a central role in regulating climate evolution.
While it is well-known that the global stratification profile exhibits a
strong dependence to Southern Ocean dynamics and in particular to wind and
buoyancy forcing, we demonstrate here that the stratification is also acutely
sensitive to the mesoscale eddy energy dissipation time-scale. Within the
context of a global ocean circulation model with an energy constrained
mesoscale eddy parameterization, it is shown that modest variations in the
eddy energy dissipation time-scale lead to significant variations in key
metrics relating to ocean circulation, namely the Antarctic Circumpolar
Current transport, Atlantic Meridional Overturning Circulation strength, and
global ocean heat content, over long time-scales. The results highlight a need
to constrain uncertainties associated with eddy energy dissipation for climate
model projections over centennial time-scales, but also for paleoclimate
simulations over millennial time-scales.

Geophysical Research Letters

Department of Physics, University of Oxford Department of Ocean Science, Hong
Kong University of Science and Technology Center for Ocean Research in Hong
Kong and Macau, Hong Kong University of Science and Technology Sorbonne
Universités (University Pierre et Marie Curie Paris 6)-CNRS-IRD-MNHN, LOCEAN
Laboratory, Paris School of Mathematics and Maxwell Institute for Mathematical
Sciences, The University of Edinburgh

Julian<EMAIL_ADDRESS>
Key metrics of global ocean circulation (ACC transport, AMOC strength, OHC
anomaly) acutely sensitive to eddy energy dissipation time-scale.

Modest variations in the dissipation time-scale has an comparable effect to
significant variations in the Southern Ocean wind forcing.

Constraints on the dissipation time-scale critical to long-time integrations
of ocean climate models such as paleoclimate scenarios.

## Plain Language Summary

Modest uncertainties in mesoscale eddy energy dissipation time-scale translate
to significant variations in the global ocean circulation and heat content on
long time-scales. A 50% change in the eddy energy dissipation time-scale has a
similar effect to halving and doubling of Southern Ocean wind stress in quasi-
equilibrium calculations, leading to a global ocean heat content change that
is an order of magnitude larger than those typically found in modern day era
reconstructions and projections. The results highlight a need to combine
theoretical, modeling, and observational efforts to constrain the
uncertainties in eddy energy dissipation for climate projections and
paleoclimate reconstructions.

##  1 Overview and key findings

Evolution of the ocean stratification plays a fundamental role in climate
evolution, through the associated consequences for the global meridional
overturning circulation. Reconstructions of past climate together with the use
of numerical models have highlighted how shoaling and weakening of the
Atlantic Meridional Overturning Circulation (AMOC), associated with changes in
the deep/abyssal stratification, have important consequences for the global
energy, oxygen and carbon cycles [<]e.g.,¿ZhangVallis13, Adkins13, Ferrari-et-
al14, Burke-et-al15, Bopp-et-al17, Jansen17, Takano-et-al18,
GalbraithdeLavergne19. In particular, the Southern Ocean is “
_disproportionately important_ ” [Newman . (2019)] for the global
stratification profiles because of the connection in the stratification
profiles, with the implication that Southern Ocean processes can exert a
control on the global overturning circulation.

It is known that the Southern Ocean stratification is primarily dependent on
wind forcing [Toggweiler Samuels (1995), Toggweiler . (2006)], buoyancy
forcing [Hogg (2010), Morrison . (2011), Jansen (2017)], and to eddy dynamics
[Munday . (2013), Farneti . (2015), Bishop . (2016)]. Focusing on mesoscale
eddies, an extra complication arises since there are notable divergences in
model response depending on how mesoscale eddies are represented, between
whether they are represented explicitly or parameterized [Munday . (2013),
Farneti . (2015), Bishop . (2016)], and the form of the parameterization
[Hofman Morales Maqueda (2011), Viebahn Eden (2012), Meredith . (2012), Munday
. (2013), Farneti . (2015), Bishop . (2016)]. While the issue of mesoscale
eddy representation “ _frequently have a larger effect on ocean climate
sensitivity than the total effect of other classes of parameterizations_ ”
[Fox-Kemper . (2019)], there have been advances on the eddy parameterization
aspect, where the role of eddy energy in mesoscale eddy parameterizations is
increasingly being studied [Eden Greatbatch (2008), Marshall Adcroft (2010),
Marshall . (2012), Eden . (2014)]. Models with parameterized eddies employing
eddy energy constrained eddy diffusivities or transport coefficients display
improved model responses that are closer to the responses displayed in
analogous high resolution models [Jansen Held (2014), Jansen, Held . (2015),
Jansen, Adcroft . (2015), Mak . (2017), Mak . (2018), Klöwer . (2018), Bachman
(2017)]. In particular, the GEOMETRIC parameterization [Marshall . (2012), Mak
. (2017), Marshall . (2017), Mak . (2018)] — effectively rescaling the
standard Gent–McWilliams [Gent McWilliams (1990), Gent . (1995)] eddy
transport coefficient by the total eddy energy according to rigorous
mathematical identities [Marshall . (2012), Maddison Marshall (2013)] and
supported in diagnoses of eddy resolving calculations [Bachman . (2017)] —
imparts an Antarctic Circumpolar Current and Atlantic Meridional Overturning
Circulation sensitivity to changes in the Southern Ocean wind forcing to
idealized ocean climate models that are closer the analogous high resolution
models [Mak . (2018)].

A natural question to ask is how strong is the influence of mesoscale eddy
dynamics on the Southern Ocean as well as global ocean circulation. Since
there is a link between eddy energy and the degree of feedback arising from
the eddies, the mesoscale eddy energy dissipation time-scale should therefore
play an important role: if more energy is drained from the mesoscale eddy
field, the associated eddy form stress weakens, thereby reducing vertical
momentum transport, in turn modifying the momentum balance. Following this
line of argument, it was argued in Marshall-et-al17 that the overall
circumpolar transport in the Southern Ocean should increase with increasing
eddy energy dissipation, and so, by thermal wind shear relation, lead to
steeply tilting isopycnals in the Southern and ocean and a deepening of the
global pycnocline depth over long time-scales.

The extent of the influence of eddy energy dissipation on the global
circulation is the primary focus of the present work. A key finding here is
that a 50% variation around a control time-scale has a similar effect to
halving and doubling of the present day Southern Ocean wind stress on the
modeled Antarctic Circumpolar Current transport, Atlantic Meridional
Overturning Circulation strength, and the global integrated ocean heat content
anomaly (Fig. 1 and Fig. 3), attributed primarily to changes in the global
pycnocline depth. While the Southern Ocean wind forcing is not expected to
vary so dramatically, the extent of plausible mesosacle eddy energy
dissipation time-scale is not known, due to a lack of theoretical and
observation constraints currently available. The results here thus highlight a
crucial need to combine theoretical, modeling, and observational efforts to
constrain the uncertainties in eddy energy dissipation, not only from a
theoretical point of view for understanding, but also for practical purposes
in constraining uncertain model parameters for numerical models used in
climate projections and paleoclimate reconstructions.

##  2 Method and model description

The principal focus here is on quasi-equilibrium sensitivities of the global
overturning circulation to the eddy energy dissipation time-scale. While one
might consider employing an eddy resolving ocean model for such a study, the
associated computational costs are prohibitive. Thus we employ a model with
parameterized eddies, and utilize the Nucleus for European Modelling of the
Ocean (NEMO, v3.7dev r8666) [Madec (2008)] in the global configuration (ORCA)
with realistic bathymetry, employing the tri-polar ORCA grid [Madec Imbard
(1996)] and the LIM3 ice model [Rousset . (2015)]. The present ORCA1 model has
a nominal horizontal resolution of $1^{\circ}$, employs 46 uneven vertical
levels, and is initialized with WOA13 climatology [Locarnini . (2013), Zweng .
(2013)]. The model employs the TEOS-10 equation of state [Roquet . (2015)],
with the atmospheric forcing modeled by the NCAR bulk formulae with normal
year forcing [Large Yeager (2009)]. Sea surface salinity but not temperature
restoration is included to reduce model drift.

An energetically constrained mesoscale eddy parameterization scheme is
required, and for our investigation the GEOMETRIC parameterization for
mesoscale eddies [Marshall . (2012), Mak . (2018)] was chosen and implemented
in NEMO (see Supplementary Information for implementation details). Briefly,
GEOMETRIC computes a horizontally and temporally varying coefficient for eddy
induced advection [Gent McWilliams (1990), Gent . (1995)] according to [<]cf.
Eq. 4 of¿Mak-et-al18

$\kappa_{\rm gm}=\alpha\frac{\int
E\;\mathrm{d}z}{\int(M^{2}/N)\;\mathrm{d}z},$ (1)

where $M$ and $N$ are the horizontal and vertical buoyancy frequencies,
$\alpha$ is a non-dimensional tuning parameter (bounded in magnitude by 1),
and $E$ is the total (potential and kinetic) eddy energy. The depth-integrated
eddy energy $\int E\;\mathrm{d}z$ is provided by a parameterized eddy energy
budget given by [<]cf. Eq. 2 of¿Mak-et-al18

$\frac{\mathrm{d}}{\mathrm{d}t}\int
E\;\mathrm{d}z+\underbrace{\nabla_{H}\cdot\left(\left(\widetilde{\mathbf{u}}^{z}-|c|\,\mathbf{e}_{x}\right)\int
E\;\mathrm{d}z\right)}_{\textnormal{advection}}=\underbrace{\int\kappa_{\rm
gm}\frac{M^{4}}{N^{2}}\;\mathrm{d}z}_{\textnormal{source}}-\underbrace{\lambda\int(E-E_{0})\;\mathrm{d}z}_{\textnormal{dissipation}}+\underbrace{\eta_{E}\nabla^{2}_{H}\int
E\;\mathrm{d}z}_{\textnormal{diffusion}}.$ (2)

The depth-integrated eddy energy is advected by the depth average flow
$\widetilde{\mathbf{u}}^{z}$ and propagated westward at the long Rossby wave
phase speed $|c|$ [Chelton . (2011), Klocker Marshall (2014)], has growth
arising from slumping of mean density surfaces, and diffused in the horizontal
[Grooms (2015), Ni, Zhai, Wang Marshall (2020), Ni, Zhai, Wang Hughes (2020)],
with $\nabla_{H}$ denoting the horizontal gradient operator and the $\eta_{E}$
the associated eddy energy diffusivity. A linear dissipation of eddy energy at
rate $\lambda$ (but maintaining a minimum eddy energy level $E_{0}$) is
utilized, so $\lambda^{-1}$ is the eddy energy dissipation time-scale. For
this work, $\alpha=0.04$ is prescribed, partially informed by the results of
Poulsen-et-al19, and $\eta_{E}=500\ \mathrm{m}^{2}\ \mathrm{s}^{-1}$ was
chosen. While the Gent–McWilliams coefficient follows the prescription given
in Eq. 1, the isoneutral diffusion coefficient [Griffies (1998)] is kept
constant at $1000\ \mathrm{m}^{2}\ \mathrm{s}^{-1}$. Values of the isopycnal
slopes used to compute the parameterized eddy energy, eddy induced advection
and isoneutral diffusion are limited to $1/100$ in the interior, and linearly
decreased from the base of the model mixed layer to zero at the surface to
maintain no flux conditions. For more details of the GEOMETRIC
parameterization and its implementation in NEMO, please see the Supporting
Information.

Given the lack of constraints on the values and uncertainties associated with
the eddy energy dissipation time-scale, for simplicity we take the eddy energy
dissipation time-scale $\lambda^{-1}$ to be a constant in space and time, with
a control value of $\lambda^{-1}=100$ days [<]e.g.,¿[Marshall & Zhai, pers.
comm.]Melet-et-al15, and consider a range of values spanning approximately 50%
about the control value (in this case six experiments with $\lambda^{-1}$
ranging 60 to 160 days in increments of 20 days). A control model utilizing
$\lambda^{-1}=100$ days was first spun up for 1500 years, after which the
perturbation experiments were integrated for a further 1600 years; see
Supplementary Information Fig. S1-4 for some of the resulting climatology. All
metrics presented in this work were diagnosed from data averaged over the last
100 model years.

##  3 Results

###  3.1 Sensitivity to eddy energy dissipation time-scale

The key metrics of interest here are the total Antarctic Circumpolar Current
(ACC) transport, Atlantic Meridional Overturning Circulation (AMOC) strength,
and the globally integrated Ocean Heat Content (OHC) anomaly relative to the
control calculation, respectively given as the transport through the model
Drake passage, the transport over the top $1000\ \mathrm{m}$ at the model 26∘
N on the Western side of the Atlantic, and the global integrated conservative
temperature multiplied accordingly by the density and heat capacity. Fig. 1
compares these metrics diagnosed from experiments varying the eddy energy
dissipation time-scale. Increasing the dissipation time-scale (i.e., decreased
damping of the eddies) leads to a substantial decrease in the ACC transport,
AMOC strength, and total ocean heat content anomaly, which can be attributed
to the deepening of the global pycnocline, consistent with theoretical
arguments [Marshall . (2017)]. In particular, we note that the changes in the
OHC anomalies found in post-industrial period reconstructions [Levitus .
(2012), Cheng . (2017), Cheng . (2019), Zanna . (2019)] are typically on the
order of $10^{23}\ \mathrm{J}$ ($100\ \mathrm{ZJ}$), while the changes to
total OHC associated with the uncertainties in eddy energy dissipation time-
scale here can be an order of magnitude larger ($10^{24}\ \mathrm{J}$).

Figure 1: Diagnostics from the varying eddy energy dissipation time-scale
experiments. Diagnostics are: ($a$) ACC transport (total in solid lines,
thermal wind component in dashed lines); ($b$) AMOC strength; ($c$) domain-
integrated ocean heat content anomalies as solid lines, where the anomalies
are relative to the control calculation (one times wind amplification and
dissipation time-scale of 100 days) with the value of 21,300 ZJ.

The distribution of the lateral depth-integrated ocean heat content for
varying eddy energy dissipation time-scale is shown in Fig. 2. Varying the
eddy energy dissipation time-scale leads to a significant global change in the
OHC anomalies, attributed mainly to the changes in pycnocline depth
(Supporting Information, Fig. S6-9). Note also that the changes appear to be
most significant over the Southern Ocean and in the Atlantic basin, attributed
to significant changes in the AMOC as well as the overturning within the
Southern Ocean (Supporting Information, Fig. S3-4).

Figure 2: Depth-integrated ocean heat content anomaly (relative to the control
calculation with dissipation time-scale of 100 days and one times wind
amplification) for varying dissipation experiments ($a,b$, at $\lambda^{-1}=$
160 and 60 days respectively).

Here, changing the eddy energy dissipation time-scale $\lambda^{-1}$ affects
the total eddy energy $E$, which in turn impacts the Gent–McWilliams
coefficient $\kappa_{\rm gm}$. While the significant changes to global OHC and
circulation arising from changing $\kappa_{\rm gm}$ has been noted before
[<]e.g.¿ZhangVallis13, the fundamental difference here is that the
sensitivities are arising through uncertainties in the eddy energy dissipation
that happens to impact the Gent–McWilliams parameter, and the eddy energy
dissipation is a process that in principle is perhaps more amenable to be
constrained by theoretical, numerical or observational means.

###  3.2 Sensitivity to Southern Ocean wind forcing

For completeness, experiments varying Southern Ocean wind forcing were also
performed. The zonal wind stress over the Southern Ocean region within the
model is amplified instead of the imposed zonal wind speed, so that any
modifications to the ocean surface evaporation and turbulent fluxes as
calculated through the bulk formulae occurs through changes to the ocean state
rather than the imposed wind forcing. Two sets of perturbation experiments
were performed: ($a$) a $\kappa_{\rm gm}$ that is varying in the horizontal
and in time as given by Eq. 1, with no further re-tuning (denoted GEOM), and
($b$) a prescribed $\kappa_{\rm gm}$ diagnosed from the last 100 years of the
control spin-up (denoted GM), but one that is now time-independent although
still spatially varying.

Fig. 3 shows the sensitivities of the same global ocean climatological metrics
to changes in the imposed Southern Ocean wind forcing. In Fig. 3($a$), while
the total ACC transport (solid lines) increases with wind forcing, the ACC
thermal wind transport (dashed lines) in the GEOM calculations (orange lines,
cross markers) is relatively insensitive to changes in the wind forcing,
demonstrating the eddy saturation phenomenon [Hallberg Gnanadesikan (2006),
Meredith Hogg (2006), Munday . (2013), Farneti . (2015), Bishop . (2016)]. On
the other hand, the corresponding AMOC strength shown in Fig. 3($b$) in the
GEOM calculations display a reduced sensitivity to changes in the Southern
Ocean wind forcing relative to the GM case, and is related to the phenomenon
of eddy compensation [Gent Danabasoglu (2011), Hofman Morales Maqueda (2011),
Viebahn Eden (2012), Meredith . (2012), Munday . (2013), Farneti . (2015),
Bishop . (2016)]. The aforementioned sensitivities coincide with a weaker
sensitivity of the global pycnocline depth to changes in Southern Ocean wind
forcing [Marshall . (2017), Mak . (2018)]. Comparing Fig. 3 with Fig. 1, a 50%
change about a chosen eddy energy dissipation time-scale has a comparable
effect to halving or doubling the Southern Ocean wind forcing in the diagnosed
metrics. As with the varying eddy energy dissipation experiments, the OHC
anomalies are particularly significant over the Southern Ocean and the
Atlantic basin (Supporting Information, Fig. S5). In the present case of
varying wind stress, however, the notable variations in the OHC anomalies are
attributed to significant changes in the abyssal watermass properties
(Supporting Information, Fig. S6-9). The observed changes in the watermass
properties may perhaps be attributed to a modified sea ice extent (Supporting
Information, Fig. S10) via changes in the sea ice export by the wind, leading
to changes in deep water formation and abyssal watermass properties, analogous
to the mechanism proposed in Ferrari-et-al14 and Burke-et-al15.

Figure 3: Diagnostics from the varying Southern Ocean wind stress experiments
(with GEOM and GM calculations in orange and blue respectively), showing:
($a$) ACC transport (total in solid lines, thermal wind component in dashed
lines); ($b$) AMOC strength; ($c$) domain-integrated ocean heat content
anomaly, where the anomalies are relative to the respective control
calculations (one times wind amplification and dissipation time-scale of 100
days) with the value of 21,300 ZJ.

##  4 Summary and outlooks

The present work demonstrates that, within the context of a global
configuration ocean model with an energetically constrained mesoscale eddy
parameterization, modest and perhaps not implausible variations in the
mesoscale eddy energy dissipation time-scale translate to significant
sensitivities of the diagnosed Antarctic Circumpolar Current transport,
Atlantic Meridional Overturning Circulation strength, and the global Ocean
Heat Content over long time-scales in the modeled ocean. The physical reasons
for the sensitivity is that modifying the eddy energy dissipation leads to
changes in the mesoscale eddy dynamics in the Southern Ocean, that in turn
lead to significant changes to the global ocean stratification over long time-
scales. The sensitivity of the aforementioned key ocean climatological metrics
to eddy energy dissipation time-scale is found to be comparable to those found
for varying Southern Ocean wind forcing, where a 50% change about a chosen
eddy energy dissipation time-scale has a comparable effect to halving or
doubling the Southern Ocean wind forcing in the diagnosed ocean circulation
metrics. In particular, changes to the globally integrated Ocean Heat Content
anomalies can vary by up to an order of magnitude larger than for
reconstructions for total ocean heat content for the anthropogenic period
[Levitus . (2012), Cheng . (2017), Cheng . (2019), Zanna . (2019)], and
comparable to the end of 21${}^{\textnormal{st}}$ projections under the
Representative Concentration Pathways scenarios [<]see figure in¿Cheng-et-
al19. While the changes in the Southern Ocean wind forcing are not expected to
vary to the extent considered in this work [<]e.g.¿Lin-et-al18, there are no
strong theoretical, numerical or observational constraints on the eddy energy
dissipation time-scale and its distribution (but see next paragraph on works
towards constraining the energy fluxes of the contributing processes). There
is thus a need to combine and dedicate theoretical, modeling and observational
efforts to constrain the uncertainties in the eddy energy dissipation time-
scale, given the impact the associated uncertainties can have.

In the present work the eddy energy dissipation is linear [<]cf.,¿Klymak18
with a time-scale that is a prescribed constant in space and time. One
particular consequence of a prescribed spatially constant eddy energy
dissipation time-scale may be seen in Fig. 4, which shows the total (kinetic
and potential) eddy energy diagnosed from a high-resolution global
configuration model and from the control experiment here. While the eddy
energy signature displays some similarities in terms of spatial patterns in
the Southern Ocean and Western Boundary Current regions, there is clearly room
for improvement for the parameterized case. For example, the eddy energy
signature in the parameterized case is too weak in the Western Boundary
Currents and in the equatorial region, attributed to the fact that the
spatially constant eddy energy dissipation time-scale was chosen somewhat with
the Southern Ocean in mind, and is probably too short for the ocean basin
regions (Marshall & Zhai, pers. comm.). The mesoscale eddy energy dissipation
time-scale is expected to be a more complicated function than the choice taken
here and, fundamentally, should depend on a wide variety of dynamical
processes such as bottom drag [<]e.g.¿Sen-et-al08, Ruan-et-al21, non-
propagating form drag [Klymak (2018), Klymak . (2021)], return to mean-flow
[<]e.g.¿Bachman19, Jansen-et-al19 scattering into internal waves
[<]e.g.¿Nikurashi-et-al13, Melet-et-al15, MacKinnon-et-al17, Yang-et-al18,
Sutherland-et-al19, loss of balance [<]e.g.¿Molemaker-et-al05, Barkan-et-al17,
Rocha-et-al18, and eddy killing by the wind [<]e.g.¿Xu-et-al16, Rai-et-al21.
Though the challenges in constraining the uncertainties in the eddy energy
dissipation time-scale are formidable, the observed/diagnosed eddy energy
signature can perhaps act as a target towards efforts to constrain the
aforementioned unknowns, highlighting the potential for further research
relating to ocean energetic pathways and its consequences for climate
evolution (see for example Ruan-et-al21 for a recent review of research
relating to ocean eddy energy pathways).

Figure 4: Depth-integrated total (kinetic and potential) eddy energy density
(in units of $\mathrm{J}\ \mathrm{m}^{-2}$), diagnosed from ($a$) the high
resolution ORCA0083-N01 calculation with explicit eddies, and ($b$) from the
GEOM control calculation.

While the present results are based on the choice of utilizing the GEOMETRIC
parameterization for mesoscale eddies [Mak . (2018)], given the link of the
eddy energy dissipation time-scale and the resulting eddy induced circulation
through the eddy energy (in this work through the Gent–McWilliams coefficient
for eddy induced advection), the sensitivities of key ocean climatological
metrics to the eddy energy dissipation is expected to carry over if other eddy
energy based parameterization schemes for mesoscale eddies [<]e.g.¿Jansen-et-
al15b, Bachman19 are utilized, or if analogous experiments are carried out in
global eddy resolving models varying possible mechanisms of mesoscale eddy
energy dissipation [<]e.g. bottom drag, cf.¿Marshall-et-al17, although the
magnitude of the sensitivities may differ. While the present work focused on
quasi-equilibrium calculations, similar conclusions but with reduced
magnitudes of the sensitivities are expected for centennial time-scale
calculations. What is clear, however, is that the present work has significant
consequences for paleoclimate simulations involving the ocean, such as
Paleoclimate Modelling Intercomparison Project calculations [<]PMIP,¿Kageyama-
et-al18, given the long time-scales inherently required for the related
simulations. Potential impact assessment for climate projections and
paleoclimate simulations in light of the present work will be investigated and
reported in due course.

###### Acknowledgements.

This work was funded by the UK Natural Environment Research Council grant
NE/R000999/1 and utilized the ARCHER UK National Supercomputing Service. JM
also acknowledges financial support from the RGC Early Career Scheme 2630020
and the Center for Ocean Research in Hong Kong and Macau, a joint research
center between the Qingdao National Laboratory for Marine Science and
Technology and Hong Kong University of Science and Technology. The data used
for generating the plots in this article is available through
`http://dx.doi.org/10.5281/zenodo.5732755`. The authors would like to thank
Andrew Coward for providing access to the ORCA0083-N01 dataset through the
ARCHER RDF service, Xiaoming Zhai for discussions relating to eddy energy
data, and George Nurser for discussions in relation to the GEOMETRIC
implementation and outlooks.

## References

  * Adkins (2013) Adkins13Adkins, JF.  2013\.  The role of deep ocean circulation in setting glacial cliamtes The role of deep ocean circulation in setting glacial cliamtes.  Paleoceanography28539–561.  10.1002/palo.20046
  * Bachman (2017) Bachman19Bachman, SD.  2017\.  The GM+E closure: A framework for coupling backscatter with the Gent and McWilliams parameterization The GM+E closure: A framework for coupling backscatter with the Gent and McWilliams parameterization.  Ocean Modell.13685–106.  10.1016/j.ocemod.2019.02.006
  * Bachman . (2017) Bachman-et-al17Bachman, SD., Marshall, DP., Maddison, JR. Mak, J.  2017\.  Evaluation of a scalar transport coefficient based on geometric constraints Evaluation of a scalar transport coefficient based on geometric constraints.  Ocean Modell.10944–54.  10.1016/j.ocemod.2016.12.004
  * Barkan . (2017) Barkan-et-al17Barkan, R., Winters, KB. McWilliams, JC.  2017\.  Stimulated imbalance and311the enhancement of eddy kinetic energy dissipation by internal waves Stimulated imbalance and311the enhancement of eddy kinetic energy dissipation by internal waves.  J. Phys. Oceanogr.47181–198.  10.1175/JPO-D-16-0117.1
  * Bishop . (2016) Bishop-et-al16Bishop, SP., Gent, PR., Bryan, FO., Thompson, AF., Long, MC. Abernathey, RP.  2016\.  Southern Ocean overturning compensation in an eddy-resolving climate simulation Southern Ocean overturning compensation in an eddy-resolving climate simulation.  J. Phys. Oceanogr.461575–1592.  10.1175/JPO-D-15-0177.1
  * Bopp . (2017) Bopp-et-al17Bopp, L., Resplandy, L., Untersee, A., Le Mezo, P. Kageyama, M.  2017\.  Ocean (de)oxygenation from the Last Glacial Maximum to the twenty-first century: insights from Earth System models Ocean (de)oxygenation from the Last Glacial Maximum to the twenty-first century: insights from Earth System models.  Phil. Trans. R. Soc. A37520160323.  10.1098/rsta.2016.0323
  * Burke . (2015) Burke-et-al15Burke, A., Stewart, AL., Adkins, JF., Ferrari, R., Jansen, MF. Thompson, AF.  2015\.  The glacial mid-depth radiocarbon bulge and its implications for the overturning circulation The glacial mid-depth radiocarbon bulge and its implications for the overturning circulation.  Paleoceanography301021–1039.  10.1002/2015PA002778
  * Chelton . (2011) Chelton-et-al11Chelton, DB., Schlax, MG. Samelson, RM.  2011\.  Global observations of nonlinear mesoscale eddies Global observations of nonlinear mesoscale eddies.  Prog. Oceanog.91167–216.  10.1016/j.pocean.2011.01.002
  * Cheng . (2019) Cheng-et-al19Cheng, L., Abraham, J., Hausfather, Z. Trenberth, KE.  2019\.  How fast are the oceans warming? How fast are the oceans warming?  Sci. Adv.363128–129.  10.1126/science.aav7619
  * Cheng . (2017) Cheng-et-al17Cheng, L., Trenberth, KE., Fasullo, J., Boyer, T., Abraham, J. Zhu, J.  2017\.  Improved estimates of ocean heat content from 1960 to 2015 Improved estimates of ocean heat content from 1960 to 2015.  Sci. Adv.3e1601545.  10.1126/sciadv.1601545
  * Eden . (2014) Eden-et-al14Eden, C., Czeschel, L. Olbers, D.  2014\.  Toward energetically consistent ocean models Toward energetically consistent ocean models.  J. Phys. Oceanogr.443160–3184.  10.1016/JPO-D-13-0260.1
  * Eden Greatbatch (2008) EdenGreatbatch08Eden, C. Greatbatch, RJ.  2008\.  Towards a mesoscale eddy closure Towards a mesoscale eddy closure.  Ocean Modell.20223–239.  10.1016/j.ocemod.2007.09.002
  * Farneti . (2015) Farneti-et-al15Farneti, R., Downes, SM., Griffies, SM., Marsland, SJ., Behrens, E., Bentsen, M.Yeager, SG.  2015\.  An assessment of Antarctic Circumpolar Current and Southern Ocean meridional overturning circulation during 1958-2007 in a suite of interannual CORE-II simulations An assessment of Antarctic Circumpolar Current and Southern Ocean meridional overturning circulation during 1958-2007 in a suite of interannual CORE-II simulations.  Ocean Modell.9484–120.  10.1016/j.ocemod.2015.07.009
  * Ferrari . (2014) Ferrari-et-al14Ferrari, R., Jansen, MF., Adkins, JF., Burke, A., Stewart, AL. Thompson, AF.  2014\.  Antarctic sea ice control on ocean circulation in present and glacial climates Antarctic sea ice control on ocean circulation in present and glacial climates.  Proc. Natl Acad. Sci. USA111248753–8758.  10.1073/pnas.1323922111
  * Fox-Kemper . (2019) FoxKemper-et-al19Fox-Kemper, B., Adcroft, AJ., Böning, CW., Chassignet, EP., Curchitser, EN., Danabasoglu, G.Yeager, SG.  2019\.  Challenges and prospects in ocean circulation models Challenges and prospects in ocean circulation models.  Front. Mar. Sci.665.  10.3389/fmars.2019.00065
  * Galbraith de Lavergne (2019) GalbraithdeLavergne19Galbraith, E. de Lavergne, C.  2019\.  Response of a comprehensive climate model to a broad range of external forcings: relevant for deep ocean ventilation and the development of late Cenozoic ice ages Response of a comprehensive climate model to a broad range of external forcings: relevant for deep ocean ventilation and the development of late Cenozoic ice ages.  Clim. Dyn.52623–679.  10.1007/s00382-018-4157-8
  * Gent Danabasoglu (2011) GentDanabasoglu11Gent, PR. Danabasoglu, G.  2011\.  Response to increasing southern hemisphere winds in CCSM4 Response to increasing southern hemisphere winds in CCSM4.  J. Climate244992–4998.  10.1175/JCLI-D-10-05011.1
  * Gent McWilliams (1990) GentMcWilliams90Gent, PR. McWilliams, JC.  1990\.  Isopycnal mixing in ocean circulation models Isopycnal mixing in ocean circulation models.  J. Phys. Oceanogr.20150–155.  10.1175/1520-0485(1990)020¡0150:IMIOCM¿2.0.CO;2
  * Gent . (1995) Gent-et-al95Gent, PR., Willebrand, J., McDougall, TJ. McWilliams, JC.  1995\.  Parameterizing eddy-induced tracer transports in ocean circulation models Parameterizing eddy-induced tracer transports in ocean circulation models.  J. Phys. Oceanogr.25463-474.  10.1175/1520-0485(1995)025¡0463:PEITTI¿2.0.CO;2
  * Griffies (1998) Griffies98Griffies, SM.  1998\.  The Gent–McWilliams skew flux The Gent–McWilliams skew flux.  J. Phys. Oceanogr.28831–841.  10.1175/1520-0485(1998)028¡0831:TGMSF¿2.0.CO;2
  * Grooms (2015) Grooms15Grooms, I.  2015\.  A computational study of turbulent kinetic energy transport in barotropic turbulence on the $f$-plane A computational study of turbulent kinetic energy transport in barotropic turbulence on the $f$-plane.  Phys. Fluids27101701.  10.1063/1.4934623
  * Hallberg Gnanadesikan (2006) HallbergGnanadesikan06Hallberg, R. Gnanadesikan, A.  2006\.  The role of eddies in determining the structure and response of the wind-driven Southern Hemisphere overturning: Results from the Modeling Eddies in the Southern Ocean (MESO) projects The role of eddies in determining the structure and response of the wind-driven Southern Hemisphere overturning: Results from the Modeling Eddies in the Southern Ocean (MESO) projects.  J. Phys. Ocenogr.362232–2252.  10.1175/JPO2980.1
  * Hofman Morales Maqueda (2011) HofmanMoralesMaqueda11Hofman, M. Morales Maqueda, MA.  2011\.  The response of Southern Ocean eddies to increased midlatitude westerlies: A non-eddy resolving model study The response of Southern Ocean eddies to increased midlatitude westerlies: A non-eddy resolving model study.  Geophys. Res. Lett.38L03605.  10.1029/2010GL045972
  * Hogg (2010) Hogg10Hogg, AM.  2010\.  An Antarctic Circumpolar Current driven by surface buoyancy forcing An Antarctic Circumpolar Current driven by surface buoyancy forcing.  Geophys. Res. Lett.37L23601.  10.1029/2010GL044777
  * Jansen (2017) Jansen17Jansen, MF.  2017\.  Glacial ocean circulation and stratification explained by reduced atmospheric temperature Glacial ocean circulation and stratification explained by reduced atmospheric temperature.  Proc. Natl Acad. Sci. USA114145–50.  10.1073/pnas.1610438113
  * Jansen . (2019) Jansen-et-al19Jansen, MF., Adcroft, A., Khani, S. Kong, H.  2019\.  Toward an energetically consistent, resolution aware parameterization of ocean mesoscale eddies Toward an energetically consistent, resolution aware parameterization of ocean mesoscale eddies.  J. Adv. Model. Earth Syst.11–17.  10.1029/2019MS001750
  * Jansen, Adcroft . (2015) Jansen-et-al15bJansen, MF., Adcroft, AJ., Hallberg, R. Held, IM.  2015\.  Parameterization of eddy fluxes based on a mesoscale energy budget Parameterization of eddy fluxes based on a mesoscale energy budget.  Ocean Modell.9228–41.  10.1016/j.ocemod.2015.05.007
  * Jansen Held (2014) JansenHeld14Jansen, MF. Held, IM.  2014\.  Parameterizing subgrid-scale eddy effects using energetically consistent backscatter Parameterizing subgrid-scale eddy effects using energetically consistent backscatter.  Ocean Modell.8036–48.  10.1016/j.ocemod.2014.06.002
  * Jansen, Held . (2015) Jansen-et-al15aJansen, MF., Held, IM., Adcroft, AJ. Hallberg, R.  2015\.  Energy budget-based backscatter in an eddy permitting primitive equation model Energy budget-based backscatter in an eddy permitting primitive equation model.  Ocean Modell.9215–26.  10.1016/j.ocemod.2015.07.015
  * Kageyama . (2018) Kageyama-et-al18Kageyama, M., Braconnot, P., Harrison, SP., Haywood, AM., Jungclaus, JH., Otto-Bliesner, BL.Zhou, T.  2018\.  The PMIP4 contribution to CMIP6 – Part 1: Overview and over-arching analysis plan The PMIP4 contribution to CMIP6 – Part 1: Overview and over-arching analysis plan.  Geosci. Model Dev.,1131033–1057.  10.5194/gmd-11-1033-2018
  * Klocker Marshall (2014) KlockerMarshall14Klocker, A. Marshall, DP.  2014\.  Advection of baroclinic eddies by depth mean flow Advection of baroclinic eddies by depth mean flow.  Geophys. Res. Lett.41L060001.  10.1002/2014GL060001
  * Klöwer . (2018) Klower-et-al18Klöwer, M., Jansen, MF., Claus, M., Greatbatch, RJ. Thomsen, S.  2018\.  Energy budget-based backscatter in a shallow water model of a double gyre basin Energy budget-based backscatter in a shallow water model of a double gyre basin.  Ocean Modell.1321–11.  10.1016/j.ocemod.2018.09.006
  * Klymak (2018) Klymak18Klymak, J.  2018\.  Non-propagating form drag and turbulence due to stratified flow over large-scale abyssal hill topography Non-propagating form drag and turbulence due to stratified flow over large-scale abyssal hill topography.  J. Phys. Oceanogr.482383–2395.  10.1175/JPO-D-17-0225.1
  * Klymak . (2021) Klymak-et-al21Klymak, J., Balwada, D., Naveira Garabato, AC. Abernathey, R.  2021\.  Parameterizing nonpropagating form drag over rough bathymetry Parameterizing nonpropagating form drag over rough bathymetry.  J. Phys. Oceanogr.511489–1501.  10.1175/JPO-D-20-0112.1
  * Large Yeager (2009) LargeYeager09Large, WG. Yeager, S.  2009\.  The global climatology of an interannually varying air-sea flux data set The global climatology of an interannually varying air-sea flux data set.  Clim. Dynam.33341–364.  10.1007/s00382-008-0441-3
  * Levitus . (2012) Levitus-et-al12Levitus, S., Antonov, JI., Boyer, TP., Baranova, OK., Garcia, HE., Locarnini, RA.Zweng, MM.  2012\.  Work ocean hat content and thermosteric sea level change (0–2000 m), 1955–2010 Work ocean hat content and thermosteric sea level change (0–2000 m), 1955–2010.  Geophys. Res. Lett.39L10603.  10.1029/2012GL051106
  * Lin . (2018) Lin-et-al18Lin, X., Zhai, X., Wang, Z. Munday, DR.  2018\.  Mean, variability, and trend of Southern Ocean wind stress: Role of wind fluctuations Mean, variability, and trend of Southern Ocean wind stress: Role of wind fluctuations.  J. Climate313557–3573.  10.1175/JCLI-D-17-0481.1
  * Locarnini . (2013) Locarini-et-al13Locarnini, RA., Mishonov, AV., Antonov, JI., Boyer, TP., Garcia, HE., Baranova, OK.Seidov, D.  2013\.  World Ocean Atlas 2013, Volume 1: Temperature World Ocean Atlas 2013, Volume 1: Temperature.  S. Levitus (), NOAA Atlas NESDIS 73 Noaa atlas nesdis 73 ( 1–40).
  * MacKinnon . (2017) MacKinnon-et-al17MacKinnon, JA., Alford, MH., Ansong, JK., Arbic, BK., Barna, A., Briegleb, BP.Zhao, Z.  2017\.  Climate process team on internal-wave driven ocean mixing Climate process team on internal-wave driven ocean mixing.  J. Fluid Mech.982429–2454.  10.1175/BAMS-D-16-0030.1
  * Maddison Marshall (2013) MaddisonMarshall13Maddison, JR. Marshall, DP.  2013\.  The Eliassen–Palm flux tensor The Eliassen–Palm flux tensor.  J. Fluid Mech.72969–102.  10.1017/jfm.2013.259
  * Madec (2008) Madec-NEMOMadec, G.  2008\.  NEMO ocean engine NEMO ocean engine.  Note du Pôle de modélisation, Institut Pierre-Simon Laplace (IPSL), No. 27.
  * Madec Imbard (1996) MadecImbard96Madec, G. Imbard, M.  1996\.  A global ocean mesh to overcome the North Pole singularity A global ocean mesh to overcome the North Pole singularity.  Clim. Dyn.12381–388.  10.1007/BF00211684
  * Mak . (2018) Mak-et-al18Mak, J., Maddison, JR., Marshall, DP. Munday, DR.  2018\.  Implementation of a geometrically informed and energetically constrained mesoscale eddy parameterization in an ocean circulation model Implementation of a geometrically informed and energetically constrained mesoscale eddy parameterization in an ocean circulation model.  J. Phys. Oceanogr.482363–2382.  10.1175/JPO-D-18-0017.1
  * Mak . (2017) Mak-et-al17aMak, J., Marshall, DP., Maddison, JR. Bachman, SD.  2017\.  Emergent eddy saturation from an energy constrained parameterisation Emergent eddy saturation from an energy constrained parameterisation.  Ocean Modell.112125–138.  10.1016/j.ocemod.2017.02.007
  * Marshall Adcroft (2010) MarshallAdcroft10Marshall, DP. Adcroft, AJ.  2010\.  Parameterization of ocean eddies: Potential vorticity mixing, energetics and Arnold’s first stability theorem Parameterization of ocean eddies: Potential vorticity mixing, energetics and Arnold’s first stability theorem.  Ocean Modell.321571–1578.  10.1016/j.ocemod.2010.02.001
  * Marshall . (2017) Marshall-et-al17Marshall, DP., Ambaum, MHP., Maddison, JR., Munday, DR. Novak, L.  2017\.  Eddy saturation and frictional control of the Antarctic Circumpolar Current Eddy saturation and frictional control of the Antarctic Circumpolar Current.  Geophys. Res. Lett.44286–292.  10.1002/2016GL071702
  * Marshall . (2012) Marshall-et-al12Marshall, DP., Maddison, JR. Berloff, PS.  2012\.  A framework for parameterizing eddy potential vorticity fluxes A framework for parameterizing eddy potential vorticity fluxes.  J. Phys. Oceanogr.42539–557.  10.1175/JPO-D-11-048.1
  * Melet . (2015) Melet-et-al15Melet, A., Hallberg, R., Adcroft, A., Nikurashin, M. Legg, S.  2015\.  Energy flux into internal lee waves: Sensitivity to future climate changes using linear theory and a climate model Energy flux into internal lee waves: Sensitivity to future climate changes using linear theory and a climate model.  J. Climate282365–2384.  10.1175/JCLI-D-14-00432.1
  * Meredith . (2012) Meredith-et-al12Meredith, MP., Garabato, ACN., Hogg, AM. Farneti, R.  2012\.  Sensitivity of the overturning circulation in the Southern Ocean to decadal changes in wind forcing Sensitivity of the overturning circulation in the Southern Ocean to decadal changes in wind forcing.  J. Climate2599–110.  10.1175/2011JCLI4204.1
  * Meredith Hogg (2006) MeredithHogg06Meredith, MP. Hogg, AM.  2006\.  Circumpolar response of Southern Ocean eddy activity to a change in the Southern Annular Mode Circumpolar response of Southern Ocean eddy activity to a change in the Southern Annular Mode.  Geophys. Res. Lett.33L16608.
  * Molemaker . (2005) Molemaker-et-al05Molemaker, MJ., McWilliams, JC. Yavneh, I.  2005\.  Baroclinic instability and loss of balance Baroclinic instability and loss of balance.  J. Phys. Oceanogr.351505–1517.  10.1175/JPO2770.1
  * Morrison . (2011) Morrison-et-al11Morrison, AK., Hogg, AM. Ward, ML.  2011\.  Sensitivity of the Southern Ocean overturning circulation to surface buoy- ancy forcing Sensitivity of the Southern Ocean overturning circulation to surface buoy- ancy forcing.  Geophys. Res. Lett.38L14602.  10.1029/2011GL048031
  * Munday . (2013) Munday-et-al13Munday, DR., Johnson, HL. Marshall, DP.  2013\.  Eddy saturation of equilibrated circumpolar currents Eddy saturation of equilibrated circumpolar currents.  J. Phys. Oceanogr.43507–532.  10.1175/JPO-D-12-095.1
  * Newman . (2019) Newman-et-al19Newman, L., Heil, P., Trebilco, R., Katsumata, K., Constable, A., van Wijk, E.Spreen, G.  2019\.  Delivering sustained, coordinated, and integrated observations of the Southern Ocean for global impact Delivering sustained, coordinated, and integrated observations of the Southern Ocean for global impact.  Front. Mar. Sci.6433.  10.3389/fmars.2019.00433
  * Ni, Zhai, Wang Hughes (2020) Ni-et-al20bNi, Q., Zhai, X., Wang, G. Hughes, CW.  2020\.  Widespread mesoscale dipoles in the global ocean Widespread mesoscale dipoles in the global ocean.  J. Geophys. Res. Oceans125e2020JC016479.  10.1029/2020JC016479
  * Ni, Zhai, Wang Marshall (2020) Ni-et-al20aNi, Q., Zhai, X., Wang, G. Marshall, DP.  2020\.  Random movement of mesoscale eddies in the global ocean Random movement of mesoscale eddies in the global ocean.  J. Phys. Oceanogr.502341–2357.  10.1175/JPO-D-19-0192.1
  * Nikurashin . (2013) Nikurashi-et-al13Nikurashin, M., Vallis, GK. Adcroft, A.  2013\.  Routes to energy dissipation for geostrophic flows in the Southern Ocean Routes to energy dissipation for geostrophic flows in the Southern Ocean.  Nature Geosci.6L08610.  10.1038/NGEO1657
  * Poulsen . (2019) Poulsen-et-al19Poulsen, MB., Jochum, M., Maddison, JR., Marshall, DP. Nuterman, R.  2019\.  A geometric interpretation of Southern Ocean eddy form stress A geometric interpretation of Southern Ocean eddy form stress.  J. Phys. Oceanogr.492553–2570.  10.1175/JPO-D-18-0220.1
  * Rai . (2021) Rai-et-al21Rai, S., Hecht, M., Maltrud, ME. Aluie, H.  2021\.  Scale of oceanic eddy killing by wind from global satellite observations Scale of oceanic eddy killing by wind from global satellite observations.  Sci. Adv.7eabf4920.  10.1126/sciadv.abf4920
  * Rocha . (2018) Rocha-et-al18Rocha, CB., Wagner, GL. Young, WR.  2018\.  Stimulated generation: extraction of energy from balanced flow by near-inertial waves Stimulated generation: extraction of energy from balanced flow by near-inertial waves.  J. Fluid Mech.847417–451.  10.1017/jfm.2018.308
  * Roquet . (2015) Roquet-et-al15aRoquet, F., Madec, G., McDougall, TJ. Barker, PM.  2015\.  Accurate polynomial expressions for the density and specific volume of seawater using the TEOS-10 standard Accurate polynomial expressions for the density and specific volume of seawater using the TEOS-10 standard.  Ocean Modell.9029–43.  10.1016/j.ocemod.2015.04.002
  * Rousset . (2015) Rousset-et-al15Rousset, C., Vancoppenolle, M., Madec, G., Fichefet, T., Flavoni, S., Barthélemy, A.Vivier, F.  2015\.  The Louvain-La-Neuve sea ice model LIM3.6: global and regional capabilities The Louvain-La-Neuve sea ice model LIM3.6: global and regional capabilities.  Geosci. Model Dev.82991–3005.  10.5194/gmd-8-2991-2015
  * Ruan . (2021) Ruan-et-al21Ruan, X., Wenegrat, JO. Gula, J.  2021\.  Slippery bottom boundary layers: The loss of energy from the general circulation by bottom drag Slippery bottom boundary layers: The loss of energy from the general circulation by bottom drag.  Geophys. Res. Lett.48e2021GL094434.  10.1029/2021GL094434
  * Sen . (2008) Sen-et-al08Sen, A., Scott, RB. Arbic, BK.  2008\.  Global energy dissipation rate of deep-ocean low-frequency flows by quadratic bottom boundary layer drag: Computations from current-meter data Global energy dissipation rate of deep-ocean low-frequency flows by quadratic bottom boundary layer drag: Computations from current-meter data.  Geophys. Res. Lett.35L09606.  10.1029/2008GL033407
  * Sutherland . (2019) Sutherland-et-al19Sutherland, BR., Achatz, U., Caulfield, CP. Klymak, JM.  2019\.  Recent progress in modeling imbalance in the atmosphere and ocean Recent progress in modeling imbalance in the atmosphere and ocean.  Phys. Rev. Fluids4010501.
  * Takano . (2018) Takano-et-al18Takano, Y., Ito, T. Deutsch, C.  2018\.  Projected centennial oxygen trends and their attribution to distinct ocean climate forcings Projected centennial oxygen trends and their attribution to distinct ocean climate forcings.  Global Biogeochem. Cycles321329–1349.  10.1029/2018GB005939
  * Toggweiler . (2006) Toggweiler-et-al06Toggweiler, JR., Russel, JL. Carson, SR.  2006\.  Midlatitude westerlies, atmospheric CO2, and climate change during the ice ages Midlatitude westerlies, atmospheric CO2, and climate change during the ice ages.  Paleoceanography21PA2005.  10.1029/2005PA001154
  * Toggweiler Samuels (1995) ToggweilerSamuels95Toggweiler, JR. Samuels, B.  1995\.  Effect of Drake passage on the global thermohaline circulation Effect of Drake passage on the global thermohaline circulation.  Deep-Sea Res., Part 1, Oceanogr. Res. Pap.424477–500.  10.1016/0967-0637(95)00012-U
  * Viebahn Eden (2012) ViebahnEden12Viebahn, J. Eden, C.  2012\.  Standing eddies in the Meridional Overturning Circulation Standing eddies in the Meridional Overturning Circulation.  J. Phys. Oceanogr.421496–1508.  10.1175/JPO-D-11-087.1
  * Xu . (2016) Xu-et-al16Xu, C., Zhai, X. Shang, XD.  2016\.  Work done by atmospheric winds on mesoscale ocean eddies Work done by atmospheric winds on mesoscale ocean eddies.  Geophys. Res. Lett.4312174–12180.  10.1002/2016GL071275
  * Yang . (2018) Yang-et-al18Yang, L., Nikurashin, M., Hogg, AM. Sloyan, BM.  2018\.  Energy loss from transient eddies due to lee wave generation in the Southern Ocean Energy loss from transient eddies due to lee wave generation in the Southern Ocean.  J. Phys. Oceanogr.482867–2885.  10.1175/JPO-D-18-0077.1
  * Zanna . (2019) Zanna-et-al19Zanna, L., Khatiwala, S., Gregory, JM., Ison, J. Heimbach, P.  2019\.  Global reconstruction of historical ocean heat storage and transport Global reconstruction of historical ocean heat storage and transport.  Proc. Natl. Acad. Sci. USA1161126–1131.  10.1073/pnas.1808838115
  * Zhang Vallis (2013) ZhangVallis13Zhang, Y. Vallis, GK.  2013\.  Ocean heat uptake in eddying and non-eddying ocean circulation models in a warming climate Ocean heat uptake in eddying and non-eddying ocean circulation models in a warming climate.  J. Phys. Oceanogr.432211–2229.  10.1175/JPO-D-12-078.1
  * Zweng . (2013) Zweng-et-al13Zweng, MM., Reagan, JR., Antonov, JI., Locarnini, RA., Mishonov, AV., Boyer, T.Biddle, MM.  2013\.  World Ocean Atlas 2013, Volume 2: Salinity World Ocean Atlas 2013, Volume 2: Salinity.  S. Levitus (), NOAA Atlas NESDIS 74 Noaa atlas nesdis 74 ( 1–39).