# Empirical Analysis of Limits for Memory Distance
in Recurrent Neural Networks

Steffen Illium1, Thore Schillman1, Robert Müller1, Thomas Gabor 1 and Claudia
Linnhoff-Popien 1
1Institute for Informatics, LMU Munich, Oettingenstr. 67, Munich, Germany
<EMAIL_ADDRESS>

###### Abstract

Common to all different kinds of recurrent neural networks (RNNs) is the
intention to model relations between data points through time. When there is
no immediate relationship between subsequent data points (like when the data
points are generated at random, e.g.), we show that RNNs are still able to
remember a few data points back into the sequence by memorizing them by heart
using standard backpropagation. However, we also show that for classical RNNs,
LSTM and GRU networks the distance of data points between recurrent calls that
can be reproduced this way is highly limited (compared to even a loose
connection between data points) and subject to various constraints imposed by
the type and size of the RNN in question. This implies the existence of a hard
limit (way below the information-theoretic one) for the distance between
related data points within which RNNs are still able to recognize said
relation.

##  1 INTRODUCTION

Neural networks (NNs) have established themselves as the premier approach to
approximate complex functions [Engelbrecht, 2007]. Although NNs introduce
comparatively many free parameters (i.e., weights that need to be assigned to
correctly instantiate the function approximation), they can still be trained
efficiently using backpropagation [Rumelhart et al., 1986], often succeeding
the performance of more traditional methods. This makes neural networks rather
fit for tasks like supervised learning and reinforcement learning, when hard-
coded algorithms are complicated to build or possess long processing times. In
general, NNs can be trained using backpropagation whenever a gradient can be
computed on the loss function (i.e., the internal function is differentiable).

One of the crucial limitations of standard NNs is that their architecture is
highly dependent on rather superficial properties of the input data. To change
the architecture during training, one has to deviate from pure gradient based
training methods [Stanley and Miikkulainen, 2002]. However, these methods
still do not allow the shape of input data to change dynamically, thus
allowing various input lengths to be processed by the same NN model.

Common use cases like the classification of (various length) time series data
have inspired the concept of _Recurrent Neural Networks_ (RNNs, cf. [Rumelhart
et al., 1986]): They apply the network to a stream of data piece by piece,
while allowing the network to connect to the state it had when processing the
previous piece. Thus, information can be carried over when processing the
whole input stream. We provide more details on the inner mechanics of RNNs and
their practical implementation in Section 2. _Long Short-Term Memory networks_
(LSTMs, cf. [Hochreiter and Schmidhuber, 1997]) are an extension of RNNs: They
introduce more connections (and thus also more free parameters) that can
handle the passing of state information from input piece to input piece more
explicitly (again, cf. Section 2 for a more detailed description). Based on
the LSTM architecture, Cho et al. introduced the _Gated Recurrent Unit_ (GRU)
with smaller numbers of free parameters but comparable training achievements
depending on the task [Cho et al., 2014, Chung et al., 2014].

All kinds of RNNs have been built around the idea that carrying over some
information throughout processing the input data stream allows the networks to
recognize connections and correlations between data points along the stream.
For example, when we want to predict the weather by the hour, it might be
beneficial to carry over data from the last couple of days as to know how the
weather was at the same time of day yesterday. However, it is also clear that
(as the number of connections ”to the past” is still fixed given a single
instance of an RNN) the amount of information that can be carried over is
subject to a hard theoretical limit. In practice, of course, NNs are hardly
capable of precisely manipulating every single bit of their connections, so
the practical limit for remembering past pieces of the input stream must be
way below the inputs vector information memory size.

In this paper, we first show that the memorization of random inputs can in
fact be learned by RNNs and then recognize a hard practical limit to how far a
given RNN architecture can remember the past. We construct the memorization
experiment in Section 3). Our experiments show that a distance limit, i.e.,
the loss of the ability to store and reproduce values, occurs very
prominently, of course varying with the type and architecture of the RNN (cf.
Section 4).

While these experiments explore the basic workings of RNNs, we are not aware
of other research assessing a comparable limitation on RNNs. Closest related
work is discussed in Section 5. The existence of a hard limit in memory
distance may have several implications on the design and usage of RNNs in
theory and practice like how big an RNN needs to be so that it is able to
remember the weather from 24 hours ago, e.g., which we briefly discuss as this
paper concludes (cf. Section 6).

##  2 FOUNDATIONS

Figure 1: ”A recurrent neural network and the unfolding in time of the
computation involved in its forward computation.” Image taken from [LeCun et
al., 2015].

Figure 2: ”Illustration of (a) LSTM and (b) gated recurrent units. (a) $i$,
$f$ and $o$ are the input, forget and output gates, respectively. $c$ and
$\tilde{c}$ denote the memory cell and the new memory cell content. (b) $r$
and $z$ are the reset and update gates, and $h$ and $\tilde{h}$ are the
activation and the candidate activation.” Image taken from [Chung et al.,
2014].

We consider the three most popular types of RNNs for experiments on their
respective memory limits. In this section, we give a quick recap on their most
relevant properties.

###  2.1 Recurrent Neural Networks

_Recurrent Neural Networks_ (RNNs) were the first models introduced to process
data sets that imply a temporal dependency on sequential inputs [Schmidhuber,
2015].

While (given an input vector $\mathbf{x}=\langle x_{i}\rangle_{1\leq i\leq
q,i\in\mathbb{N}}$ of length $q\in\mathbb{N}$) one element $x_{i}$ at a time
is processed, a hidden state vector $\mathbf{h}$ carries the history of all
the past elements of the sequence [LeCun et al., 2015]. Finally, a tanh
activation function shifts the internal states to the nodes of the output
layer. Given a sequential input $\mathbf{x}$ and an activation function $g$,
the output state for every time step $t$ is $h_{t}=g(Wx_{t}+Uh_{t-1})$ [Chung
et al., 2014]. RNNs can be thought of as single layer of neural cells, which
is reused multiple times on a sequence while taking track of all subsequent
computations (cp. Figure 2). Therefore, every position of the output vector
$z$ is influenced by the combination of every previously evaluated input as
described by:

$z=tanh((x\odot wx+b)+(ws\odot s))$ (1)

###  2.2 Long Short-Term Memory Networks

While regular RNNs are known for their problems in understanding long-time
dependencies as the gradient either explodes or vanishes [LeCun et al., 2015,
Chung et al., 2014, Schmidhuber, 2015], Hochreiter and Schmidhuber introduced
the  _Long Short-Term Memory (LSTM)_ unit [Hochreiter and Schmidhuber, 1997].
The main difference lies in the ”connection to itself at the next time step
through a memory cell” [LeCun et al., 2015] as well as a forget gate, that
regulates the temporal connection (cp. Figure 2, a). Layers equipped with LSTM
cells can learn what to forget while remembering valuable inputs in long
temporal sequences. However, this functionality comes at a high computational
cost in total trainable parameters $N=4\cdot(mn+n^{2}+n)$ where $m$ is the
input dimension, and $n$ is the cell dimension. This constitutes a four-fold
increase in comparison to a RNN [Lu and Salem, 2017].

###  2.3 Gated Recurrent Unit Networks

Cho et al. proposed a new LSTM variant, Gated Recurrent Unit (GRU) with
reduced parameters and no memory cell. At first ”the reset gate $r$ decides
whether the previous hidden state is ignored”, then ”the update gate $z$
selects whether the hidden state is to be updated”. In combination with the
exposure of the hidden state, the internal decision gates could be reduced
from four to ”only two gating units”, so that the overall number of trainable
parameters per layer activation decreases (see [Cho et al., 2014, Chung et
al., 2014]; cp. Figure 2, a). In short, a GRU cell is a hidden layer unit that
”includes a reset gate and an update gate that adaptively control how much
each hidden unit remembers or forgets while reading/generating a sequence”
[Cho et al., 2014]. Even though [Chung et al., 2014] did not prefer GRU over
LSTM in an empirical testing setup, they saw that the usage ”may depend
heavily on the dataset and corresponding task”.

###  2.4 Echo State Networks

_Echo state networks_ (ESNs) are a type of _reservoir computer_ ,
independently developed by Herbert Jaeger [Jaeger, 2001] and Wolfgang Maass
[Maass et al., 2002]. Reservoir computing is a framework which generalises
many neural network architectures utilising recurrent units. ESNs differ from
”vanilla” RNNs in that the neurons are sparsely connected (typically around 1%
connectivity) and the reservoir weights are fixed. Also, connections from the
input neurons to the output as well as looping connections in the reservoir
(hidden layer) are allowed, see Figure 3.

Figure 3: An ESN with all possible kinds of connections. The legend denotes
which ones are fixed, trainable and optional.

Unlike most recurrent architectures in use today, ESNs are not trained via
_backpropagation though time_ (BPTT), but simple linear regression on the
weights to the output units. Having fixed reservoir weights makes the
initialization even more important. Well initialized ESNs display the _echo
state property_ , which guarantees, after an initial transient, that the
activation state is a function of the input history:

$h_{t}=E(...,x_{t-1},x_{t})$ (2)

Here, $E$ denotes an echo function $E=(e_{1},...,e_{N})$ with
$e_{i}:U^{\mathbb{-N}}\rightarrow\mathbb{R}$ being the echo function of the
$i$-th unit.

Several conditions for echo states have been proposed, the most well-known
stating that the spectral radius and norm of the reservoir matrix must be
below unity [Jaeger, 2001]. This has been shown to neither be sufficient nor
necessary [Yildiz et al., 2012, Mayer, 2017], but remains a good rule of
thumb.

The amount of information from past inputs that can be stored in the ESNs
transient state dynamics is often called its short-term memory and is
typically measured by the _memory capacity_ (MC), proposed by Jaeger [Jaeger,
2001]. It is defined by the coefficient of correlation between input and
output, summed over all time delays

$MC=\sum_{k=1}^{\infty}MC_{k}=\sum_{k=1}^{\infty}\frac{cov^{2}(x(t-k),y_{k}(t))}{Var(x(t))*Var(y_{k}(t))}$
(3)

and thus, measures how much of the inputs variance can be recovered by
optimally trained output units. It has been shown by Jaeger that the MC of an
ESN with linear output units is bounded above by its number of reservoir
neurons N.

ESNs predate the deep learning revolution. Nonetheless, deep variants of ESNs
have recently been proposed [Gallicchio and Micheli, 2017] to make use of the
advantages of layering. Reservoir computing is an active field of research, as
the models are less expensive to train, often biologically plausible
[Buonomano and Merzenich, 1995] and allow for precise description of state
dynamic properties due to the fixed weights of the reservoir.

##  3 METHODOLOGY

To reveal the limits of distance in which information can be carried over from
past states $h_{t-n}$ through recurrent connections of RNN type networks, we
construct a task in which the RNN is forced to remember a specific input of
the input data stream $\mathbf{x}$.

This implies that there must not be any correlation between the pieces of the
input data stream so that there is no way to reconstruct past pieces without
fully passing in their values through the recurrent connections to all future
applications of the RNN. We implemented this experiment using random-length
sequences of random numbers as input data to the RNN. The task (cf. Def. 1) is
the reproduction of the $p$th-from-last input throughout a training setup for
a set $p\in\mathbb{N}$.

However, since the sequences are of random lengths, the RNN cannot recognize
which piece of the input stream will be needed until the very last
application. This forces the RNN to carry over all the information it can from
the input stream. We argue that this is the canonical experiment to test
memorization abilities of RNNs. This task definitions are in stark contrast to
common applications of RNNs where one explicitly expects some correlations
between the pieces of the input sequence and thus uses RNNs to approximate a
model of said connections.

######  Definition 1 (Random Memorization Task).

Let $q$ be a random number uniformly sampled from
$[10;15]\subseteq\mathbb{N}$. Let $\mathbf{x}=\langle{x_{i}}\rangle_{1\leq
i\leq q,i\in\mathbb{N}}$ be a sequence of random numbers $x_{i}$ each
uniformly sampled from $[0;1]\subseteq\mathbb{R}$. Let $\mathcal{R}$ be an RNN
that is applied to the sequence $\mathbf{x}$ in order with the result
$\mathcal{R}(\mathbf{x})$ of the last application. The _memorization task of
the $p$th-from-last position_ is given by the minimizing goal function
$f_{p}(\mathcal{R},\mathbf{x})=|x_{q+1-p}-\mathcal{R}(\mathbf{x})|$.

Note that the random length of the sequence is substantial. Otherwise, an RNN
trained to reproduce the $p$th-from-last position may count the position of
the pieces of the input sequence from the start, making it rather easy to
remember the $q$th-from-last (i.e., first) position. In complexity terms, this
task requires space complexity of O(1). These claims are supported by our
experiments in Section 4, where we also test a fixed-length memorization task:

######  Definition 2 (Fixed-Length Random Memorization Task).

Let $q=10$. For all else, refer to Def. 1.

##  4 EXPERIMENTAL RESULTS

To evaluate the task described in Definitions 1 we conducted several
experiments with different configurations of RNNs, LSTMs, and GRUs. A
configuration is specified by the number of layers $l$ of the network and the
number of cells per layer $c$. We disregarded any complex architectural
pattern by just examining standard fully connected recurrent layers as well as
modern network training techniques such as gradient clipping or weight decay
in favor of comparability. Layers’ outputs from more than $c\geq 1$ cells are
added up. Weight updates are applied by _stochastic gradient descent_ (SGD)
and the truncated backpropagation algorithm. Training data is generated by
sampling from a Mersenne Twister pseudo-random number generator. Every network
configuration was trained for up to 5000 epochs for 5 different random seeds
if not stated otherwise.

Random Memorization We start by evaluating the random memorization task (cf.
Definition 1). We run a simple scenario with $l=1$ layer and $c=5$ cells
within that layer. We conducted $10$ independent training runs each for every
position to memorize from ”position -1”, i.e., the last position or $p=1$, to
”position -10”, i.e., the $10$th-from-last position or $p=10$. Figure 4 shows
the results for an RNN plotting the loss over training time: Reproducing
positions $1\leq p\leq 4$ is learned quickly so that most runs overlap at the
bottom of the plot. $p=5$ seems to be harder to learn as some runs take
distinctively more time to reach a loss close to zero, but the $5$th-from-last
random memorization task is still learned reliably. Positions $6\leq p\leq 8$
seem to stagnate not being able to reduce the loss beyond around $0.15$, which
still means that some information about these positions’ original value is
being passed on, though. However, the network greatly loses accuracy and
cannot reliably reproduce the values. We observe another case in most runs of
positions $p=9$ and $p=10$, where the loss seems to stagnate at around $0.25$:
Note that for the random memorization task of sequence values
$x_{i}\in[0;1]\subseteq\mathbb{R}$ a static network without any memorization
can simply always predict the fixed value $\mathcal{R}(\mathbf{x})=0.5$ to
achieve an average loss of $0.25$ on a sequence of random numbers. At least,
that trick is easily learned by all setups of RNNs.

We compare those results to experiments executed for LSTMs in Figure 5. First,
we notice that the LSTM takes much longer training time for the simpler
positions $1\leq p\leq 5$, which are all learned eventually but show a clear
separation. This result is intuitive considering the additional number of
weights (total size $4\cdot(mn+n^{2}+n)$ for $m$ input elements and $n$
neurons), i.e., trainable parameters, for the same size of network compared to
RNNs (total size $mn+n^{2}+n$). On the other hand, we can observe that LSTMs
still show some learning progress on all positions $p$ after the 500 epochs
plotted here, especially having a run of position $p=6$ taking off towards
near-optimal loss. Thus, there is reason to believe that LSTMs may eventually
prove to be considerably more powerful than RNNs given significantly longer
training times on our task.

Figure 4: Absolute loss for an RNN with $l=1$ layer and $c=5$ cells when asked
to reproduce the $p$th-from-last position (for $p\in[1;10]$) of a random
sequence of random length $q\in[10;15]$. Every color shows $10$ different
training runs (with different random seed) for a single $p$. Figure 5:
Absolute loss for an LSTM with $l=1$ layer and $c=5$ cells when asked to
reproduce the $p$th-from-last position (for $p\in[1;10]$) of a random sequence
of random length $q\in[10;15]$. Every color shows $10$ different training runs
(with different random seed) for a single $p$.

For a more complete overview of the limits of performance in the random
memorization task we ran this experiment for different number of layers
$l\in[1;5]$ and cells per layer $c\in[1;20]\subseteq\mathbb{N}$ as well as for
the $p$th-from-last positions $p\in[1;20]\subseteq\mathbb{N}$ as a target.
Figure 6 shows the results for RNNs. The color of the respective dots encoding
the average loss over $5$ independent runs, the scale of the loss is again
fitted to $0.25$ as all RNNs manage to learn that fixed guess as discussed
above. We can observe that the $1$-layer $1$-cell network at the bottom left
only learns to reproduce the last value of the random sequence, which it can
do without any recurrence by just learning the identity function. It is thus
clear why all RNNs manage to solve this instance. Furthermore, we can clearly
see a monotone increase in difficulty when increasing the position, i.e.,
forcing the network to remember more of the sequence’s past. This shows that
no network can learn a generalization of the task that can memorize arbitrary
positions as presumed when constructing the task.

What is possibly the most interesting result of these experiments is that the
number of instances resulting in average losses of around $[0.05;0.2]$
(colored in light blue, green) is relatively small: Certain positions can
either be reproduced very clearly or barely better than the fixed guess
(purple, orange), cf. Figure 7. This implies a strong connection between the
amount of information that can be passed on throughout the RNN. Furthermore,
we can observe that more cells and more layers (which also bring more cells to
the table) clearly help in solving the memorization task. Interestingly, even
parallel architectures as shown for the networks with only one layer ($l=1$)
become better with more cells. The connection here appears almost linear:
Within our result data, a one-layer network with $c=n$ cells can always
memorize the $n$th-from-last position ($p=n$) very well.

Figure 6: Mean absolute loss for various configurations of RNNs with
$l\in[1;5]$ layers with $c\in[1;20]$ cells each when asked to reproduce the
$p$th-from-last position (for $p\in[1;20]$) of a random sequence of random
length $q\in[10;15]$. Every data point is averaged over $5$ runs. Figure 7:
Standard Deviation for various configurations of RNNs with $l\in[1;5]$ layers
with $c\in[1;20]$ cells each when asked to reproduce the $p$th-from-last
position (for $p\in[1;20]$) of a random sequence of random length
$q\in[10;15]$. Every data point is averaged over $5$ runs.

Figure 8 shows the same experimental setup as described above for LSTM
networks. In line with the observation about the loss plots (cf. Figure 5),
LSTMs do not perform as well as RNNs on this task, i.e., not satisfyingly
solving any position $p=10$ instance ($\textit{mean\\_absolute\\_loss}\geq
0.04$). However, as described above we conjecture this stems from the fact
that the higher number of weights present in an LSTM network requires
substantially longer training times beyond what could be performed as part of
this experiment. These results still have practical implications: As this
relatively small task consumes large computational resources, LSTMs seem to be
inefficient at least at (smaller instances of) the random memorization task,
at least for single layer networks.

We also conducted the same experiments for GRU networks, as shown in Figure 8.
The GRUs’ behavior appears most similar to LSTMs’. We will use this experiment
as a basis to refrain from showing GRU plots for other experiments since GRUs
have performed that similarly throughout all our experiments. However, note
that GRUs require fewer parameters ($N=3\cdot(mn+n^{2}+n)$), making them
interesting for efficient practical applications [Dey and Salemt, 2017].

Figure 8: Mean absolute loss for various configurations of LSTMs (left) and
GRU (right) with $l\in\\{1,2,3\\}$ layers with $c\in[1;10]$ cells each when
asked to reproduce the $p$th-from-last position (for $p\in[1;10]$) of a random
sequence of random length $q\in[10;15]$. Every data point is averaged over
$10$ runs.

Our results indicate that Jaegers upper bound on ESN short-term memory remains
valid in the context of deep RNN architectures trained by backpropagation
(RNN, LSTM, GRU). For stacked, multi-layer architectures, the maximum memory
capacity could not be reached. Instead, we observed an upper bound of $MC\leq
N-(l-1)$ with $N=c*l$.

###  4.1 Fixed-Length Random Memorization

We now quickly discuss the scenario of Definition 2: In this case, we fix the
length of the random sequence $q=10$. This means that in theory it is possible
to ”count” the positions from the start and thus know exactly when the
currently processed piece of the input data stream needs to be memorized,
sparing the effort of memorizing as many of them as possible until the end. We
thus consider this task strictly easier and our experimental results shown in
Figure 9 support this. However, the strategy of ”counting from the start”
seems to be subject to some memory limitations as well. Thus, while every
instance that is easily learned in the case of variable-length sequences is
also learnable in this scenario, further learnable instances appear from the
right of the plot, i.e., $p=10$ for the 1-layer 2-cells network and $p\geq 8$
for the 2-layers 2-cells-each network. This shows that the RNN is not exactly
able to count through the sequence it is given but able to memorize pieces
from the start, within a similar limit as from the end of the sequence. As
LSTMs and GRUs show a similar behavior, we omit the respective plots for
brevity.

Figure 9: Mean absolute loss for various configurations of RNNs with
$l\in\\{1,2,3\\}$ layers with $c\in[1;10]$ cells each when asked to reproduce
the $p$th-from-last position (for $p\in[1;10]$) of a random sequence of fixed
length $q=10$. Every data point is averaged over $10$ runs.

##  5 RELATED WORK

Recent studies on recurrent neural networks mostly focused on memory afforded
by attractor dynamics or weight changes while learning. Both are only useful
in the context of temporally correlated data, as the model settling in a low
dimensional state space would not be desirable when every data point is
equally important. As RNNs have been conceived precisely to recognize and
model dependencies in sequential data, their ability to deal with i.i.d. or
uncorrelated input might seem to be of little importance. We argue that this
capability is an integral part of RNN memory and performance, which has long
been recognized in the reservoir computing literature. For tasks where the
output should depend strongly on the whole input history, we require a unique,
clearly distinguishable hidden state for every input. Using uncorrelated data,
we can approximate the total amount of information the architecture can store
in its transient activation dynamics.

In addition, and specific to RNN-architectures a lot of empirical evaluations
on various LSTM and RNN configurations and tasks already exist. Hochreiter and
Schmidhuber for example evaluated LSTM cells from the beginning. They
conjectured that gradient-based approaches suffer from practical inability to
precisely count discrete time steps and therefore assume the need for an
additional counting mechanism [Hochreiter and Schmidhuber, 1997]. However, our
experiments show that an explicit counting mechanism might not be necessary,
even for models, trained on random sequences (i.e., pure noise). At least for
the extend of our observations.

First experiments on the ability of an RNN to carry information over some
distance in the input sequence were conducted [Cleeremans et al., 1989].
Likening an RNN to a finite state automaton (for more details on that line of
thought [Giles et al., 1992]), they trained a network to recognise an
arbitrary but fixed regular grammar. Interestingly, they concluded that a
recurrent network can encode information about long-distance contingencies
only as long as information about critical past events is relevant at each
time step, which again contrasts our results on random sequences.

Little research has been performed on random sequences of variable length as
inputs for RNNs. One could argue that experiments 4 and 5 by Hochreiter and
Schmidhuber are of relevance here [Hochreiter and Schmidhuber, 1997]. However,
we like to stress that our work is focused on random sequences both in length
and value to reproduce an element $x_{i}$ on a specific position without
providing any form of hint such as additional labelling or encoding. We see
our work in contrast to these early experiments, as we test for maximal signal
remoteness for a given computational budget and not necessary for memory
endurance in general. In other words, given an input vector of randomly
generated length, they ask for the binary answer to whether a specific symbol
came before another. This tests if the gradient can uphold for enough
recurrent steps so that the two bits of information can be recognized. In
contrast, given a random input vector of randomly generated length, we ask the
networks to reproduce the $p$th-last entry. We also considered classic RNNs
that do not possess gated memory vectors.

Most recently Gonon et al. analysed the memory and forecasting capacities of
linear and nonlinear RNN (for independent and non-independent inputs). They
stated an upper bound for both while generalizing Jaegers statements regarding
MC estimates. [Gonon et al., 2020]

Additionally, Bengio et al. performed related experiments to investigate the
problem of vanishing and exploding gradients. Furthermore, their findings
explain how robust the transport of information is under the influence of
noisy and unrelated information. However, they did not consider pure
memorization ”by heart”. [Bengio et al., 1994]

As one might confuse the benchmark problems in this paper with work on the
topic of generalization and over-fitting, we would like to point out, that our
targeted notion of capacity is very different. Although unrelated in
methodology (they used regular fully connected feed-forward networks), Zhang
et al. provide an analysis quite similar in wording: They consider a network’s
capacity as the amount of input-output pairs from the training set that can be
exactly reproduced by the network (after training) [Zhang et al., 2016]. This
kind capacity is measured on a whole data set of input vectors. We consider
memory capacity within a single call of a trained network for a single input
vector that is passed in piece by piece, which is an entirely different setup.

Neural Turing Machines currently are a promising extension of classical NNs
designed specifically to offer more powerful operations to manage a hidden
state ([Graves et al., 2014]). We reckon that the Random Memorization Task is
trivial to them. The same applies for neural networks architectures that
involve attention mechanism, like the very prominent Transformer [Vaswani et
al., 2017].

##  6 CONCLUSION

We defined the Random Memorization Task for RNNs. Despite it contradicting the
intuition behind the usage of recurrent neural networks, classical RNNs, LSTM
and GRU networks were able to memorize a random input sequence to a certain
extent, depending on configuration and architecture. There is a discernible
borderline between past positions in the sequence that could be memorized and
those that could not, which correlates with the MC formulated by Jaeger. We
therefor conclude, that Jaegers MC formula is applicable for calculating the
memory limit in respect to the RNN’s type and architecture.

While our experiments are very limited in scale (due to the already
computationally expensive nature of the experiments), we observe the trend
that more cells increase the memory limit. The limiting factor that we
observed, at least for vanilla-RNN’s, is the VEGP. However, it is important to
note that for current RNN applications the ratio of input sequence length and
cells inside the RNN is usually in favor of the number of cells, so that we
would not expect the memory limit to play an important role in the application
of the typical RNN. Still, we hope this research represents an additional step
towards a theoretical framework concerned with the learnability of problems
using specific machine learning techniques.

## REFERENCES

  * Bengio et al., 1994 Bengio, Y., Simard, P., Frasconi, P., et al. (1994).  Learning long-term dependencies with gradient descent is difficult.  IEEE transactions on neural networks, 5(2):157–166.
  * Buonomano and Merzenich, 1995 Buonomano, D. V. and Merzenich, M. M. (1995).  Temporal information transformed into a spatial code by a neural network with realistic properties.  Science, 267(5200):1028–1030.
  * Cho et al., 2014 Cho, K., Van Merriënboer, B., Gulcehre, C., Bahdanau, D., Bougares, F., Schwenk, H., and Bengio, Y. (2014).  Learning phrase representations using RNN encoder-decoder for statistical machine translation.  arXiv preprint arXiv:1406.1078.
  * Chung et al., 2014 Chung, J., Gulcehre, C., Cho, K., and Bengio, Y. (2014).  Empirical evaluation of gated recurrent neural networks on sequence modeling.  arXiv preprint arXiv:1412.3555.
  * Cleeremans et al., 1989 Cleeremans, A., Servan-Schreiber, D., and McClelland, J. L. (1989).  Finite state automata and simple recurrent networks.  Neural computation, 1(3):372–381.
  * Dey and Salemt, 2017 Dey, R. and Salemt, F. M. (2017).  Gate-variants of gated recurrent unit (GRU) neural networks.  In 2017 IEEE 60th International Midwest Symposium on Circuits and Systems (MWSCAS), pages 1597–1600. IEEE.
  * Engelbrecht, 2007 Engelbrecht, A. P. (2007).  Computational intelligence: an introduction.  John Wiley & Sons.
  * Gallicchio and Micheli, 2017 Gallicchio, C. and Micheli, A. (2017).  Deep echo state network (deepesn): A brief survey.  arXiv preprint arXiv:1712.04323.
  * Giles et al., 1992 Giles, C. L., Miller, C. B., Chen, D., Chen, H.-H., Sun, G.-Z., and Lee, Y.-C. (1992).  Learning and extracting finite state automata with second-order recurrent neural networks.  Neural Computation, 4(3):393–405.
  * Gonon et al., 2020 Gonon, L., Grigoryeva, L., and Ortega, J.-P. (2020).  Memory and forecasting capacities of nonlinear recurrent networks.  Physica D: Nonlinear Phenomena, 414:132721.
  * Graves et al., 2014 Graves, A., Wayne, G., and Danihelka, I. (2014).  Neural turing machines.  arXiv preprint arXiv:1410.5401.
  * Hochreiter and Schmidhuber, 1997 Hochreiter, S. and Schmidhuber, J. (1997).  Long short-term memory.  Neural Computation, 9(8):1735–1780.
  * Jaeger, 2001 Jaeger, H. (2001).  The “echo state” approach to analysing and training recurrent neural networks-with an erratum note.  Bonn, Germany: German National Research Center for Information Technology GMD Technical Report, 148(34):13.
  * LeCun et al., 2015 LeCun, Y., Bengio, Y., and Hinton, G. (2015).  Deep learning.  Nature, 521(7553):436.
  * Lu and Salem, 2017 Lu, Y. and Salem, F. M. (2017).  Simplified gating in long short-term memory (LSTM) recurrent neural networks.  In 2017 IEEE 60th International Midwest Symposium on Circuits and Systems (MWSCAS), pages 1601–1604. IEEE.
  * Maass et al., 2002 Maass, W., Natschläger, T., and Markram, H. (2002).  Real-time computing without stable states: A new framework for neural computation based on perturbations.  Neural computation, 14(11):2531–2560.
  * Mayer, 2017 Mayer, N. M. (2017).  Echo state condition at the critical point.  Entropy, 19(1):3.
  * Rumelhart et al., 1986 Rumelhart, D. E., Hinton, G. E., and Williams, R. J. (1986).  Learning representations by back-propagating errors.  nature, 323(6088):533.
  * Schmidhuber, 2015 Schmidhuber, J. (2015).  Deep learning in neural networks: An overview.  Neural networks, 61:85–117.
  * Stanley and Miikkulainen, 2002 Stanley, K. O. and Miikkulainen, R. (2002).  Evolving neural networks through augmenting topologies.  Evolutionary computation, 10(2):99–127.
  * Vaswani et al., 2017 Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, L., and Polosukhin, I. (2017).  Attention is all you need.  In NIPS.
  * Yildiz et al., 2012 Yildiz, I. B., Jaeger, H., and Kiebel, S. J. (2012).  Re-visiting the echo state property.  Neural networks, 35:1–9.
  * Zhang et al., 2016 Zhang, C., Bengio, S., Hardt, M., Recht, B., and Vinyals, O. (2016).  Understanding deep learning requires rethinking generalization.  arXiv preprint arXiv:1611.03530.


#  Venus boundary layer dynamics:
eolian transport and convective vortex

Maxence Lefèvre  Department of Physics (Atmospheric, Oceanic and Planetary
Physics), University of Oxford, Oxford, UK

(Accepted in Icarus)

## Abstract

Few spacecraft have studied the dynamics of Venus’ deep atmosphere, which is
needed to understand the interactions between the surface and atmosphere.
Recent global simulations suggest a strong effect of the diurnal cycle of
surface winds on the depth of the planetary boundary layer. We propose to use
a turbulent-resolving model to characterize the Venus boundary layer and the
impact of surface winds for the first time. Simulations were performed in the
low plain and high terrain at the Equator and noon and midnight. A strong
diurnal cycle is resolved in the high terrain, with a convective layer
reaching 7 km above the local surface and vertical wind of 1.3 m/s. The
boundary layer depth in the low plain is consistent with the observed
wavelength of the dune fields. At noon, the resolved surface wind field for
both locations is strong enough to lift dust particles and engender micro-
dunes. Convective vortices are resolved for the first time on Venus.

##  1 Introduction

The interaction between the surface and the atmosphere is a major aspect of
exchanges of heat and angular momentum, impacting the thermal and wind shear
profiles and, therefore, the atmospheric dynamics and the rotation of the
solid body itself. On Venus, however, the first 10 km above the surface remain
a largely unknown region due to the technological difficulty in probing this
region below the cloud layer.

Only a limited number of probes have been able to collect data on Venus. The
VeGa-2 probe has successfully measured the only temperature profile in that
region (Linkin et al.,, 1986). At the surface, only Venera 9 and 10 directly
measured the wind for respectively 49 min and 90 s (Avduevskii et al.,, 1977),
and several other probes like Venera 13 and 14 measured indirectly the wind
speed (Ksanfomaliti et al.,, 1983). The amplitudes of the measured wind speeds
are less than 2 m s-1 below 100 m (Lorenz,, 2016), with a higher probability
for values below 0.5 m s-1. The height of the planetary boundary layer (PBL)
and the diurnal cycle is not known, nor are the effect of the topography.
Dunes have been observed in radar measurements with Magellan (Greeley et al.,,
1992), although the lack of knowledge about the spatial and temporal
distribution of winds complicates the interpretation of dust transport.

The Institut Pierre Simon Laplace (ISPL) Venus General Circulation Model (GCM)
simulations showed the diurnal cycle of the PBL activity is correlated with
the diurnal cycle of surface winds (Lebonnois et al.,, 2018). They observed
downward katabatic winds at night and upward anabatic winds during the day
along the slopes of high-elevation terrains, resulting in a deeper PBL depth
at noon.

Yamamoto, (2011) performed turbulent-resolving simulations of the Venus PBL
with a resolution of 100 m without radiative processes and a variety of near-
surface theoretical thermal structures. This experiment yielded a PBL depth
below 2.5 km. Morellina and Bellan, (2022) studied the turbulent chemical-
species mixing at the surface, where high-density-gradient magnitude regions
are formed with larger gradients due to the supercritical conditions.

The PBL of the Earth and Mars have been studied extensively in both
observations and modelling studies. On Earth, the turbulent flux is in the
energy budget of the convective layer. The opposite is true on Mars. For
Venus, this budget is poorly known. On Earth, the PBL dynamic is associated
with the presence of the water cycle and latent heat release. The
understanding of such dynamics is crucial for the clouds and water cycle and
the energy balance of the surface (Garratt,, 1994). On Mars and Venus, the
latent heat is negligible. The depth of the Martian surface convective layer
is larger than on Earth, with also a greater temperature diurnal cycle (Hinson
et al.,, 2008). The surface slope winds have a strong impact on the thermal
structure of the PBL (Spiga et al.,, 2011).

In this study, we use Large Eddy Simulation (LES) models developed for the
clouds convective regions and to simulate the PBL convective activity at two
different locations on the surface. This is done to examine the influence of
the topography on the boundary layer convection, and two local times, noon and
midnight, to quantify the diurnal cycle.

For the first time, the radiative processes are taken into account in the
study of the Venus PBL with prescribed solar and IR heating rates. An
additional rate representing the effect of the general circulation,
heating/cooling due to large-scale wind advection, is also prescribed.

In Section 2, the model is described. The spatial and local time variabilities
of the PBL are discussed in Section 3. The PBL of Venus is compared to the
Earth and Mars in Section 4. The impact of PBL characteristics on eolian
transport is discussed in Section 5. Our conclusions are summarized in Section
6.

##  2 Modelling

###  2.1 Dynamical core

The LMD LES model is based on the dynamical core of the Advanced Research
Weather-Weather Research and Forecast (hereinafter referred to as WRF)
terrestrial model (Skamarock and Klemp,, 2008). The WRF dynamical core
integrates the fully compressible non-hydrostatic Navier-Stokes equations over
a specified area of the planet. The conservation of the mass, momentum, and
entropy are ensured by an explicitly conservative flux-form formulation of the
fundamental equations (Skamarock and Klemp,, 2008), based on mass-coupled
atmospheric variables (winds and potential temperature). The parametrization
of the unresolved small-scale eddies is carried out by a subgrid-scale
prognostic Turbulent Kinetic Energy closure by Deardorff, (1972). This
methodology has been used for extensively for Earth convection study (Moeng et
al.,, 2007), and the Martian atmosphere (Spiga et al.,, 2010), Venus cloud
convective layer (Lefèvre et al.,, 2017, 2018), and terrestrial exoplanet
convection (Lefèvre et al.,, 2021).

###  2.2 Model Physics

Due to the constant heat capacity of the dynamical core and the architecture
of the coupled radiative transfer, as well as computational time, the
radiative forcing is handled in the same way that in Lefèvre et al., (2017).
The solar and radiative heating rates are extracted from ISPL Venus GCM
(Garate-Lopez and Lebonnois,, 2018), which uses IR transfer (Lebonnois et
al.,, 2015) based on Eymet et al., (2009), with the latitudinally-varying
cloud model of Haus et al., (2014, 2015). An additional heating rate is
prescribed, representing the large-scale heating from the dynamics. In this
study, we focus on the first 10 km and the large-scale heating will come
mainly from the anabatic/katabatic slope flows. No surface and sub-surface
physics are considered in this study.

###  2.3 Simulation settings

Lebonnois et al., (2018) showed with GCM modelling that the convective depth
in the PBL was impacted by the diurnal cycle of the surface wind, and was
maximal in the steepest slope of the equatorial topographic features.
Therefore, we choose two locations at the surface with two distinct elevations
and slope environments at the Equator. Here, the incoming solar flux is
maximized in order to study the activity of the PBL where it is supposed to be
the most active. One of the locations is in the low plain, with an elevation
of -320 m at 0∘ longitude, and the other is in the western part of Ovda Regio
with an elevation of 1030 m at 80∘ longitude. The point in the plain will be
hereinafter referred to as low plain, and the point in Ovda Regio will be
referred to as high terrain. The domain of the LES simulations is flat. Due to
computational constraints, simulations of an entire Venus day were not
possible, and two local times are considered in this study: noon and midnight.
The surface is heat flux is set at 90 W m-2 at noon and -1 W m-2 at midnight
for the two locations (Lebonnois et al.,, 2018). This flux is constant in time
during the simulations over the entire domain, there is no feedback of the PBL
turbulence on the sensible flux. For the two locations, the horizontal
resolution and timestep are set at 50 m and 0.4 s. However, the size of the
surface area varies depending on local time and location, 30$\times$30 km for
the high terrain case at noon and 20$\times$20 km for the rest. The vertical
resolution and extent also depend on the location and local time, from 10 km
above the local surface with a mean resolution of 90 m for the High terrain
case at noon to 5 km above the local surface, with a resolution of 60 m. The
different horizontal domains size were chosen to allow several connective
cells in each horizontal direction, and were determined by trial and error.
The different vertical domains size were chosen to allow several kilometers
above the convective layer and were based on the Venus IPSL GCM results
(Lebonnois et al.,, 2018). To avoid the spurious reflection of gravity waves
propagating upward on the top of the model, a Rayleigh damping layer is
applied over the last 500 m with a damping coefficient of 0.01 s-1. The heat
capacity is set to a constant value over the whole domain of 1181 J K-1, a
reference value from the Venus International Reference Atmosphere (Seiff et
al.,, 1985).

Fig 1 shows the initial profiles of the temperature, potential temperature,
and the different heating rates. The temperature is colder for the high
terrain cases, and the diurnal cycle of the temperature is below 3 K for the
two location cases. At noon, a neutral layer corresponding to the convective
layer is visible below 2 km above the local surface for the low plain case,
and below 8 km above the local surface for the high terrain case. At midnight,
there is no visible neutral layer, meaning that the convective activity is
weak. Regarding the heating rates, the short wave heating is slightly greater
for the low plain case, but the thermal cooling is slightly stronger at
midnight. However, there is a strong difference for the large-scale heating
whereas for the high terrain case it is positive up to 6 km above the local
surface, with stronger values at noon. While for the low plain case, the
large-scale heating is negative in the first 1.5 km at noon and 2 km at
midnight, and then alternating between positive and negative values above.
This variability of the large-scale heating reflects the effect of the
topography and the diurnal cycle of the surface wind. This variability
reflects on the total heating rate, alternating between negative and positive
values that will enforce the convective depth. The model is initialized with
thermal profiles, winds and radiative rates that reached GCM equilibrium,
using hypotheses from the subgrid-scale parametrization that will impact the
equilibrium state of the region. With the very poor knowledge of the first 12
km, it is difficult to assess the realism of the initial state used for the
present studies, although the surface temperature and winds are consistent
with measurements. The time step of the simulations is set to 0.4 s, and the
radiative transfer is called every 100 dynamical steps, ensuring that there is
no issue with numerical precision (Rafkin and Soto,, 2020). The outputs shown
in the following sections are obtained after at least 1 Earth day of
simulations for the midnight cases, where a steady-state is reached. For the
noon cases, after 2 Earth days, there is a small temperature drift due to the
set-up of the surface flux and radiative rates, around 10-8 K/s near the
surface. A constant forcing corresponding to a single time of day is sometimes
used to study self-aggregation in Earth’s tropics (Daleu et al.,, 2015; Wing
et al.,, 2017), as well as tidally-locked rocky exoplanets (Zhang et al.,,
2017; Sergeev et al.,, 2020; Lefèvre et al.,, 2021). When no aggregation is
present, the model reaches equilibrium in a couple of days (Wing and Emanuel
2014). The Venus atmosphere can be considered dry. In the present study, there
is no water vapor and clouds, and therefore no feedback between moisture and
radiation/surface temperature. No aggregation of convection is expected. There
is no noticeable change in pressure over time in the domain for all the cases
considered. The surface wind’s amplitude range does not vary in time. The
pseudo-equilibrium is realistic enough to provide qualitative insight into the
Venus PBL dynamics.

Figure 1: Upper left : initial temperature vertical profiles (K). Upper right
: Potential temperature vertical profiles (K). Bottom left: IR (LW), solar
(SW) and large-scale (LS) heating rates vertical profiles (K/Venus day).
Bottom right : total heating rate vertical profiles (K/Venus day).

The model configuration is summarized in Table 1.

Parameter | Value
---|---
Gravity (m s-2) | 8.87
Heat Capacity (J K-1) | 1181
Surface heat flux (W m-2) | 90 (noon), -1 (midnight)
Horizontal resolution dx (m) | 50
Time step (s) | .4
Cases | Grid |  $\Delta$x and $\Delta$y (km) |  $\Delta$z (km) | dz (m)
low plain noon | 401$\times$401$\times$71 | 20 | 5.5 | 80
low plain midnight | 401$\times$401$\times$51 | 20 | 3 | 60
high terrain noon | 601$\times$601$\times$121 | 30 | 10 | 70
high terrain midnight | 401$\times$401$\times$101 | 20 | 5.5 | 55
Table 1: Planetary and atmospheric parameters for the different simulations.
$\Delta$x represents the size of the domain over the axis $x$ and dx the
resolution over the axis $x$

##  3 Spatial and temporal variability of the PBL

Fig 2 shows snapshots of the vertical and horizontal cross-sections of the
vertical wind for the high terrain at noon and midnight, and the low plain at
noon. The diurnal cycle for the high terrain location is striking: at noon the
elevation of the PBL can reach 7 km above the local surface with vertical wind
speed as high as 1.3 m s-1. Whereas at midnight, the convective layer barely
reaches 0.5 km with vertical wind speed $<$ 0.2 m s-1. This diurnal cycle is
consistent with GCM simulations (Lebonnois et al.,, 2018). The difference in
depth leads to a difference in convective cell diameter. At noon, the typical
cell diameter is around 5 km, with some cells reaching 10 km, and about 2 km
at midnight. The cellular features are elongated in the y-direction at
midnight. The Richardson number is stronger at midnight than at noon, meaning
that the shear becomes relatively strong compared to buoyancy. Such features
depend on the vertical shear of the horizontal wind, where there is very
little data to constrain this result. With similar surface flux, the elevation
of the PBL depth is different at noon, 2 km above the local surface for the
low plain compared to 7 km for the high terrain. This difference is due to the
large-scale forcing: in the high terrain this forcing is positive up to 6 km
above the local surface (Fig 1), bringing energy to the environment and
leading to a deeper convective layer. Whereas in the low plain, the large-
scale forcing is strongly negative in the first 1.5 km thus stabilizing the
atmosphere and leading to a smaller PBL depth. In the high terrain at night,
the large-scale forcing will also bring more heat, (although less than at
noon) to the environment, but the IR cooling is stabilizing the atmosphere. In
the low plain at midnight (not shown here), both the large-scale forcing and
IR cooling stabilize the atmosphere, resulting in no convective activity. This
PBL spatial variability is also consistent with the GCM simulations (Lebonnois
et al.,, 2018).

Figure 2: Snapshots of the vertical cross-section (left column row) and
horizontal cross (right column) at 100 m above the local surface of the
vertical wind (m/s) for the high terrain case at noon (top row) and at
midnight (middle row) and the low plain case at noon (bottom column).

To demonstrate the impact of the large-scale forcing on the convective depth,
simulations were carried out without large-scale forcing and wind shear, but
with the same thermal profiles. Results are shown in Fig 3. Without the large-
scale forcing, the height of the PBL decreases by several kilometers in the
high terrain and by several hundred meters in the low plain. At noon, the
large-scale circulation will heat the atmosphere. Such an impact on the PBL
depth due to the topography is not observed on Earth, but is on Mars (Spiga et
al.,, 2010). The convection-resolving model of Yamamoto, (2011) displays a
realistic surface heat flux with a PBL depth of around 2 km, consistent with
the low plain case. With higher surface flux, though unrealistic, the
convective layer can go as high as 6 km.

Figure 3: Vertical profile of the domain-average potential temperature for the
high terrain case (left) and the low plain case (right) at noon with the
large-scale heating rates and wind shear (red lines) and without the large-
scale heating rates or wind shear (blue lines).

In 3D GCMs or 1D models, the turbulence processes are not resolved, therefore
a parametrization is needed and a vertical eddy diffusivity coefficient, or
$K_{zz}$, is often used to represent convection. In LES models, like the one
presented here, the resolution is small enough to resolve the larger eddies,
and no vertical eddy diffusivity is therefore needed.

$K_{zz}$ has been estimated in the Venus atmosphere modelling, but is still
not well-constrained. The 1D models of Matsuda and Matsuno, (1978) and Takagi
et al., (2010) have estimated the vertical eddy diffusion between 102 and 103
m2 s-1. In the IPSL Venus GCM, the subgrid parametrization exhibits values as
high as 1 m2 s-1 in the high terrain (Lebonnois et al.,, 2018). Yamamoto,
(2011) convection-resolving model displays values between 1 and 104 m2 s-1.

Using Prandtl mixing-length theory (Lindzen,, 1971) the vertical eddy
diffusion can be estimated from the resolved convective plume and, as
expected, it is higher with a stronger convective activity. It reaches values
as high as 103 m2 s-1 in the high terrain at noon, consistent with the
estimation from Yamamoto, (2011) with similar PBL depth. However, the values
are several orders of magnitude higher than the IPSL Venus GCM. The
parametrization of the GCM, based on Mellor and Yamada, (1982), does not
display realistic values for vertical eddy diffusion compared to the
convective-resolving model of Yamamoto, (2011) and the present study.

##  4 Comparison with Earth and Mars

The depth of the convective layer in the low plain at noon is comparable to
the Earth’s shallow convection depth (Garratt,, 1994) but with much weaker
wind, 0.2 m s-1 for Venus, and several meters per second for Earth (Park et
al.,, 2018). The depth of the convective layer in the high terrain is somewhat
comparable to the Martin PBL depth, reaching 10 km at maximum (Hinson et al.,,
2008). This has a much lower vertical wind amplitude, around 1 m s-1 for
Venus, and 10 m s-1 for Mars (Spiga et al.,, 2010). The difference in vertical
wind is also visible in the convective velocity scale, defined as

$W^{\star}\leavevmode\nobreak\ =\leavevmode\nobreak\
\left[g.z_{h}\frac{\langle
w^{\prime}\theta^{\prime}\rangle_{max}}{\langle\theta\rangle}\right]^{1/3}$
(1)

with zh the height of the PBL. For Venus, W∗ is around 0.4 m s-1, around 2 m
s-1 for Earth (Stull,, 1988), and between 4 and 6.5 m s-1 for Mars (Spiga et
al.,, 2010). Further comparison between the three atmospheres is shown in Fig
4 with the vertical eddy heat flux and vertical velocity variance in the PBL
in dimensionless defined respectively as

$\frac{\langle w^{\prime}\theta^{\prime}\rangle_{max}}{\langle
w^{\prime}\theta^{\prime}\rangle}\hskip 28.45274pt;\hskip
28.45274pt\frac{w^{{}^{\prime}2}}{W^{\star}}$ (2)

The heat flux on Earth is maximum at the surface, whereas for Venus and Mars
fluxes are at the surface, and decrease with height. The turbulence acts to
heat the atmosphere. On Mars and Venus, the turbulent flux first increases
with height, acting to cool the atmosphere, and then decreases with height, as
on Earth, to heat the atmosphere. The low vertical wind amplitude in the PBL
shows that with a denser atmosphere, around 50 times higher than the Earth,
and 3000 times that of Mars atmospheric density, the heat transport by the
atmosphere is more efficient than on other terrestrial planets. For Venus, as
on Earth, the energy budget of the convective layer is dominated by the
convective flux, with negligible solar and IR heating rates. In the Martian
atmosphere, the radiative heating rates are dominant, with surface
temperatures that do not depend on altitude, and a pressure effect that is
obtained for similar solar heating and IR cooling rates (Spiga et al.,, 2010).
The empirical similarities in this dimensionless approach between the three
planets provide a framework for boundary layer parametrization.

Figure 4: Vertical profiles of the vertical eddy heat flux (left) and vertical
velocity variance (right) for Venus, Earth and Mars PBL. Heat flux, variance
and height are normalized respectively by the maximum vertical eddy heat flux,
the convective velocity scale W∗ and the boundary layer depth. For Venus, the
high terrain case at noon is presented, Earth data is from Stull, (1988), and
Mars data is from Spiga et al., (2010).

##  5 Eolian Transport

Only two prominent dune fields have been identified on Venus with the Magellan
radar (Greeley et al.,, 1992), Algaonice at 25∘S, 340∘E cover some 1300 km 2
with bright dunes of 0.5-5 km in length with a wavelength around 0.5 km, and
Fortuna-Meshkenet at 67∘N, 91∘E with transverse dunes of 0.5-10 km long,
0.2-0.5 km wide and spaced by an average of 0.5 km. Several sites with micro-
dunes or small-scale wind streaks have also been proposed (Weitz et al.,,
1994; Bondarenko et al.,, 2006) but the resolution of the Magellan radar is
not high enough to confirm them.

The wavelength of the observed dune fields is about 500 m. The height of the
PBL in Algaonice dune field would be slightly lower than the one at the
Equator (Fig 2-bottom line) with an averaged depth over a Venus day around 500
m, about the same order of magnitude than the wavelength of the dunes.
However, in the Fortuna-Meshkenet the averaged PBL depth would be much lower
due to the solar flux more than ten times lower, although the proximity of
Ishtar Terra could engender wind turbulence. On Earth, the wavelength of giant
dunes is controlled by the depth of the PBL (Andreotti et al.,, 2009). On
Venus, From the two Magellan dune fields observations, the relationship
between the wavelength of a giant dune and the depth of the PBL seems to be
more complicated, and more observations are needed.

The wind profile is generally approximated by the Prandlt-von Kármán equation
defining the friction velocity u⋆ as:

$u^{\star}\leavevmode\nobreak\ =\leavevmode\nobreak\
k\frac{u_{z}}{ln(z/z_{0})}$ (3)

with the Von Kármán constant, k, equals to 0.4, the surface roughness, z0, set
to 1 cm. uz is the horizontal wind at the altitude z.

Theoretical calculations and laboratory experiments estimated the threshold
friction velocity for which the dust is lifted in Venus surface condition to
be minimum around 2.5 10-2 m s-1 (Iversen et al.,, 1976), about 10 times lower
than for Earth, and 100 times lower than for Mars due to the very dense
atmosphere, depending on dust radius. Such a value depends also on the density
of the dust, data that is not well-known yet. Laboratory experiments tested
the dust transport in Venus surface conditions and showed that for surface
velocities between 0.63 and 1.5 m s-1 there was the formation of micro-dunes,
with a wavelength between 8 and 27 cm (Greeley et al.,, 1984). Below 0.63 m
s-1, the amplitude of the wind is too small to transport dust, and above 1.5 m
s-1 no dunes are formed because dust grains can be transported across the lee
of a wave onto another and therefore blurring the separation between dunes.

For both the high terrain and low plain, both at noon, the friction velocity
of the vast majority of the spatial area is above the saltation threshold,
meaning that the turbulence is able to transport dust particles, as shown in
Fig 5. On the other hand, at night for both locations, the friction velocity
is below the saltation threshold over the entire area, and no dust is
transported. At night, the horizontal wind amplitude distribution at 10 m
above the local surface for both locations is fully consistent with in-situ
measurements (Lorenz,, 2016), with wind speed $<$ 0.5 m s-1. Whereas at noon,
the majority of the wind distribution is $>$ 0.5 m s-1, with values of 1 m s-1
for the high terrain. These values are on the outer edge of the observed
values, without taking site-specific factors like slopes or local times into
account. With a deeper convective layer in the high terrain, the distribution
is broadened compared to the low plain location. In the low plain, the
horizontal wind amplitude is almost always in the amplitude range where micro-
dunes are forming, whereas in the high terrain a significant part of the
horizontal wind is too weak to transport wind to form micro-dunes. The
formation of micro-dunes is therefore more probable in the low plain at noon,
however the spatial distribution of particle reservoir is not known.

Figure 5: Snapshots of the friction velocity u⋆ (top row) and horizontal wind
amplitude distribution (bottom row) in the high terrain (left column) and low
plain (right column) at 10 m above the local surface, both at noon.

Fig 6 shows maps of the surface temperature and surface pressure in the high
terrain at noon. The temperature anomaly induced by the convection is very
small, around 0.05 K for the strongest convective activity. However, in the
pressure field, there is a significant drop, characteristic of convective
vortices at the convergence of convective cells (Rafkin et al.,, 2016).

Figure 6: Snapshots of the surface temperature (left), the surface pressure
perturbation (right) for the high terrain case at noon. The black line in the
right panel represents the cross-section in the following figure.

At these convergences, a pressure deficit is formed from rising plumes, where
the updrafts cause a vertical stretching of horizontal wind, leading to a
spiralling flow of rising air (Sinclair,, 1973; Rennó et al.,, 1998; Toigo et
al.,, 2003; Ito et al.,, 2013). The vorticity inside the vortices, visible
with the vertical component of the rotation vector $\zeta$ = ($\partial
v$/$\partial x$) - ($\partial u$/$\partial y$) with $u$ and $v$ the two
horizontal wind fields. There are vortices with both positive and negative
vortices, meaning, respectively, cyclonic and anticyclonic advection. These
patterns are also visible in the horizontal divergence $D$ = ($\partial
u$/$\partial x$) + ($\partial v$/$\partial y$), with convergence (positive
divergence values) for the cyclonic advection and divergence (negative
divergence values) for anticyclonic advection. The vortices are also regions
where the shear, $\partial u$/$\partial y$ \+ $\partial v$/$\partial x$, and
stretching, $\partial u$/$\partial x$ \- $\partial v$/$\partial y$,
deformations are strong. Fig 7 shows horizontal cross-sections of the surface
pressure, vorticity, divergence, shear and stretching at y = 25.2 km (black
line in Fig 6) where there is the presence of a strong vortex. The main vortex
is visible at x = 10 km, with the presence of two smaller vortices at x = 12.2
km and x = 12.9 km with lower pressure drop. The vortices are strongly visible
in the vorticity field, with values up to 2$\times$10-2 s-1, at least an order
of magnitude lower than Earth vortices (Spiga et al.,, 2016).The vortices are
barely discernible in the divergence fields, but generate a lot of shear and
stretching. There are generations of several strong shear and stretching
increase or drops, without noticeable pressure decreases at x = 14.7 km, 22 km
or 23.9 km, for example, due to the presence of thin convective updrafts
between convective cells.

The main vortex extends up to 5 km above the local surface, much higher than
on Earth, where it is capped by the height of the PBL (around 2 km). The two
smaller vortices only extend up to 1 km above the local surface. The lifetime
of such vortices can be several hours. No such structures are resolved in the
low plain at noon with a 50 m resolution. The convective activity is either
too weak to form vortices, or the resolution of the simulations is too coarse
to resolve thin vortices. At night for both locations, no vortices are
resolved because the convective activity is too weak to create a vortex. The
pressure drops are around 100 Pa for the largest vortices, comparable to the
highest pressure drops of dust devils on Earth (Lorenz and Lanagan,, 2014).
However, the diameters are different, almost 1 km on Venus against a few
hundreds of meters, maximum, for Earth dust-devils (Lorenz,, 2011). Convective
vortices are also measured in the Martian atmosphere (Baker et al.,, 2021),
with pressure drops being inferior to a few pascals (Nishizawa et al.,, 2016;
Spiga et al.,, 2021). The 50 m horizontal resolution used for these
simulations is sufficient to qualitatively resolve the characteristic
structure of convective vortices, but a horizontal resolution one order of
magnitude lower would be necessary to resolve a larger part of the spectrum of
kinetic energy (Spiga et al.,, 2016), as well as smaller vortices.

Figure 7: Horizontal cross-sections at y=25.2 km (black line if Fig 6) of the
surface pressure (A), the vertical component of the rotation vector (B), the
horizontal divergence of the horizontal wind (C), the shear (D) and stretching
(E) at 10 m above the local surface for the high terrain case at noon.

Lorenz, (2021) used an equation to estimate the maximum vertical wind
amplitude in a convective vortex for the different terrestrial planets and
Titan. For Venus, this maximum vertical wind amplitude was estimated at 2 m
s-1, stronger than the one in the present study around 1.3 m s-1. However, the
value of the heat capacity and surface sensitive heat flux chosen in Lorenz,
(2021) were not representative of the Venus surface at the Equator. With a
heat capacity of 1181 J K-1 and a surface sensitive heat flux of 90 W m-2, the
estimation with the equation from Lorenz, (2021) is 1.35 m s-1, very close to
the value obtained in the present model in the high terrain. The location of
the vortices is where the friction velocity is the highest and superior to the
saltation threshold. Therefore, if there is dust material to be lifted, it
will generate dust devils.

##  6 Conclusion

This study presents turbulent-resolving simulations of the Venus PBL, with for
the first a realistic thermal profile and with the impact of the large-scale
circulation. Simulations were performed at two locations on the Equator, in
the low plain and high terrain, and at two local times, noon and midnight. The
convective activity is stronger at noon due to the incoming stellar radiation.
It is also stronger in the high terrain due to the large-scale heating from
the slope winds. The vertical eddy diffusion estimated from turbulent-
resolving simulations is not consistent with the GCM, with values several
orders of magnitude higher, consistent with previous turbulent-resolving
modelling.
The PBL of Venus was compared to its equivalent in the atmosphere of the Earth
and Mars. On Venus and Mars, the turbulent flux first increases with height,
acting to cool the atmosphere, whereas on Earth the turbulence acts to heat
the atmosphere. The energy budget of the convective layer is dominated by the
convective flux, with negligible solar and IR heating rates for Venus and the
Mars. In the Martian atmosphere, the radiative heating rates are dominant.
At noon, the turbulent activity is strong enough for the friction velocity to
be above the saltation threshold needed for eolian transport. However, there
are a lot of unknowns about the composition, density, radius, and reservoir of
such particles that could drastically affect their transport. These
simulations were carried out at noon on the Equator, where the convective
activity is more intense. More simulations need to be performed at higher
latitude and other local times.
In the high terrain at noon, the strong activity in the PBL convective
vortices are resolved for the first time. Pressure drops as strong as 100 Pa
are observed, comparable to the highest pressure drops of dust devils on
Earth, but with a larger diameter, almost 1 km on Venus versus a few hundreds
km for the Earth. The vertical velocity inside such vortices can exceed 1 m
s-1. With the presence of dust, these vortices could lead to dust devils.

One of the main improvements of the model setup would be to have an
interactive surface heat flux and radiative transfer like in the Venus cloud
layer (Lefèvre et al.,, 2018), this would allow feedback between the surface
and the atmosphere.

In the IPSL Venus GCM, a nitrogen vertical gradient in the first 7 km from the
surface has been tested (Lebonnois and Schubert,, 2017; Lebonnois et al.,,
2018) and has a substantial impact on the stability of the atmosphere in this
region and leads to weaker PBL depth. Such gradient should be tested with a
similar setup in future studies.

## Acknowledgements

The author would like to thank Aymeric Spiga, Sébastien Lebonnois and Kevin
Olsen for helpful comments and discussions. The author would like to thank
Scot Rafkin and an anonymous reviewer for their help to improve the
manuscript. The author acknowledges funding from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation program
(grant agreement No. 740963/EXOCONDENSE). The author would like to acknowledge
the use of the University of Oxford Advanced Research Computing (ARC) facility
in carrying out this work. http://dx.doi.org/10.5281/zenodo.22558.

## References

  * Andreotti et al., (2009) Andreotti, B., Fourrière, A., Ould-Kaddour, F., Murray, B., and Claudin, P. (2009).  Giant aeolian dune size determined by the average depth of the atmospheric boundary layer.  Nature, 457(7233):1120–1123.
  * Avduevskii et al., (1977) Avduevskii, V. S., Vishnevetskii, S. L., Golov, I. A., Karpeiskii, I. I., Lavrov, A. D., Likhushin, V. I., Marov, M. I., Melnikov, D. A., Pomogin, N. I., and Pronina, N. N. (1977).  Measurement of wind velocity on the surface of Venus during the operation of stations Venera 9 and Venera 10.  Cosmic Research, 14(5):710–713.
  * Baker et al., (2021) Baker, M., Newman, C., Charalambous, C., Golombek, M., Spiga, A., Banfield, D., Lemmon, M., Banks, M., Lorenz, R., Garvin, J., Grant, J., Lewis, K., Ansan, V., Warner, N., Weitz, C., Wilson, S., and Rodriguez, S. (2021).  Vortex Dominated Aeolian Activity at InSight’s Landing Site, Part 2: Local Meteorology, Transport Dynamics, and Model Analysis.  Journal of Geophysical Research (Planets), 126(4):e06514.
  * Bondarenko et al., (2006) Bondarenko, N. V., Kreslavsky, M. A., and Head, J. W. (2006).  North-south roughness anisotropy on Venus from the Magellan Radar Altimeter: Correlation with geology.  Journal of Geophysical Research (Planets), 111(E6):E06S12.
  * Daleu et al., (2015) Daleu, C. L., Plant, R. S., Woolnough, S. J., Sessions, S., Herman, M. J., Sobel, A., Wang, S., Kim, D., Cheng, A., Bellon, G., Peyrille, P., Ferry, F., Siebesma, P., and van Ulft, L. (2015).  Intercomparison of methods of coupling between convection and large-scale circulation: 1. Comparison over uniform surface conditions.  Journal of Advances in Modeling Earth Systems, 7(4):1576–1601.
  * Deardorff, (1972) Deardorff, J. W. (1972).  Numerical Investigation of Neutral and Unstable Planetary Boundary Layers.  Journal of Atmospheric Sciences, 29:91–115.
  * Eymet et al., (2009) Eymet, V., Fournier, R., Dufresne, J.-L., Lebonnois, S., Hourdin, F., and Bullock, M. A. (2009).  Net exchange parameterization of thermal infrared radiative transfer in Venus’ atmosphere.  J. of Geophys. Res. (Planets), 114:E11008.
  * Garate-Lopez and Lebonnois, (2018) Garate-Lopez, I. and Lebonnois, S. (2018).  Latitudinal variation of clouds’ structure responsible for Venus’ cold collar.  Icarus, 314:1–11.
  * Garratt, (1994) Garratt, J. R. (1994).  Review: the atmospheric boundary layer.  Earth Science Reviews, 37(1):89–134.
  * Greeley et al., (1992) Greeley, R., Arvidson, R. E., Elachi, C., Geringer, M. A., Plaut, J. J., Saunders, R. S., Schubert, G., Stofan, E. R., Thouvenot, E. J. P., Wall, S. D., and Weitz, C. M. (1992).  Aeolian features on Venus - Preliminary Magellan results.  Journal of Geophysical Research, 97:13.
  * Greeley et al., (1984) Greeley, R., Marshall, J. R., and Leach, R. N. (1984).  Microdunes and other aeolian bedforms on Venus: Wind tunnel simulations.  Icarus, 60(1):152–160.
  * Haus et al., (2014) Haus, R., Kappel, D., and Arnold, G. (2014).  Atmospheric thermal structure and cloud features in the southern hemisphere of Venus as retrieved from VIRTIS/VEX radiation measurements.  Icarus, 232:232–248.
  * Haus et al., (2015) Haus, R., Kappel, D., and Arnold, G. (2015).  Radiative heating and cooling in the middle and lower atmosphere of Venus and responses to atmospheric and spectroscopic parameter variations.  Planetary and Space Science, 117:262–294.
  * Hinson et al., (2008) Hinson, D. P., Pätzold, M., Tellmann, S., Häusler, B., and Tyler, G. L. (2008).  The depth of the convective boundary layer on Mars.  Icarus, 198(1):57–66.
  * Ito et al., (2013) Ito, J., Niino, H., and Nakanishi, M. (2013).  Formation Mechanism of Dust Devil-Like Vortices in Idealized Convective Mixed Layers.  Journal of Atmospheric Sciences, 70(4):1173–1186.
  * Iversen et al., (1976) Iversen, J. D., Greeley, R., and Pollack, J. B. (1976).  Windblown dust on Earth, Mars and Venus.  Journal of Atmospheric Sciences, 33:2425–2429.
  * Ksanfomaliti et al., (1983) Ksanfomaliti, L. V., Goroshkova, N. V., and Khondyrev, V. K. (1983).  Wind velocity on the Venus surface from acoustic measurements.  Kosmicheskie Issledovaniia, 21:218–224.
  * Lebonnois et al., (2015) Lebonnois, S., Eymet, V., Lee, C., and Vatant d’Ollone, J. (2015).  Analysis of the radiative budget of the Venusian atmosphere based on infrared Net Exchange Rate formalism.  J. of Geophys. Res. (Planets), 120:1186–1200.
  * Lebonnois and Schubert, (2017) Lebonnois, S. and Schubert, G. (2017).  The deep atmosphere of Venus and the possible role of density-driven separation of CO2 and N2.  Nature Geoscience, pages 473–477.
  * Lebonnois et al., (2018) Lebonnois, S., Schubert, G., Forget, F., and Spiga, A. (2018).  Planetary boundary layer and slope winds on Venus.  Icarus, 314:149–158.
  * Lefèvre et al., (2018) Lefèvre, M., Lebonnois, S., and Spiga, A. (2018).  Three-Dimensional Turbulence-Resolving Modeling of the Venusian Cloud Layer and Induced Gravity Waves: Inclusion of Complete Radiative Transfer and Wind Shear.  Journal of Geophysical Research (Planets), 123:2773–2789.
  * Lefèvre et al., (2017) Lefèvre, M., Spiga, A., and Lebonnois, S. (2017).  Three-dimensional turbulence-resolving modeling of the Venusian cloud layer and induced gravity waves.  Journal of Geophysical Research (Planets), 122:134–149.
  * Lefèvre et al., (2021) Lefèvre, M., Turbet, M., and Pierrehumbert, R. (2021).  3D Convection-resolving Model of Temperate, Tidally Locked Exoplanets.  The Astrophysical Journal, 913(2):101.
  * Lindzen, (1971) Lindzen, R. S. (1971).  Tides and Gravity Waves in the Upper Atmosphere.  In Fiocco, G., editor, Mesospheric Models and Related Experiments, volume 25 of Astrophysics and Space Science Library, page 122\.
  * Linkin et al., (1986) Linkin, V. M., Kerzhanovich, V. V., Lipatov, A. N., Pichkadze, K. M., Shurupov, A. A., Terterashvili, A. V., Ingersoll, A. P., Crisp, D., Grossman, A. W., Young, R. E., Seiff, A., Ragent, B., Blamont, J. E., Elson, L. S., and Preston, R. A. (1986).  VEGA balloon dynamics and vertical winds in the Venus middle cloud region.  Science, 231:1417–1419.
  * Lorenz, (2011) Lorenz, R. (2011).  On the statistical distribution of dust devil diameters.  Icarus, 215(1):381–390.
  * Lorenz, (2016) Lorenz, R. D. (2016).  Surface winds on Venus: Probability distribution from in-situ measurements.  Icarus, 264:311–315.
  * Lorenz, (2021) Lorenz, R. D. (2021).  Dust devil winds: Assessing dry convective vortex intensity limits at planetary surfaces.  Icarus, 354:114062.
  * Lorenz and Lanagan, (2014) Lorenz, R. D. and Lanagan, P. D. (2014).  A Barometric Survey of Dust-Devil Vortices on a Desert Playa.  Boundary-Layer Meteorology, 153(3):555–568.
  * Matsuda and Matsuno, (1978) Matsuda, Y. and Matsuno, T. (1978).  Radiative-convective equilibrium of the Venusian atmosphere.  Meteorological Society of Japan, 56:1–18.
  * Mellor and Yamada, (1982) Mellor, G. L. and Yamada, T. (1982).  Development of a turbulence closure model for geophysical fluid problems.  Reviews of Geophysics and Space Physics, 20:851–875.
  * Moeng et al., (2007) Moeng, C., Dudhia, J., Klemp, J., and Sullivan, P. (2007).  Examining Two-Way Grid Nesting for Large Eddy Simulation of the PBL Using the WRF Model.  Monthly Weather Review, 135(6):2295–2311.
  * Morellina and Bellan, (2022) Morellina, S. and Bellan, J. (2022).  Turbulent chemical-species mixing in the Venus lower atmosphere at different altitudes: a direct numerical simulation study relevant to understanding species spatial distribution.  Icarus, 371:114686.
  * Nishizawa et al., (2016) Nishizawa, S., Odaka, M., Takahashi, Y. O., Sugiyama, K.-i., Nakajima, K., Ishiwatari, M., Takehiro, S.-i., Yashiro, H., Sato, Y., Tomita, H., and Hayashi, Y.-Y. (2016).  Martian dust devil statistics from high-resolution large-eddy simulations.  Geophysical Research Letters, 43(9):4180–4188.
  * Park et al., (2018) Park, S., Kim, S.-W., Park, M.-S., and Song, C.-K. (2018).  Measurement of Planetary Boundary Layer Winds with Scanning Doppler Lidar.  Remote Sensing, 10(8):1261.
  * Rafkin et al., (2016) Rafkin, S., Jemmett-Smith, B., Fenton, L., Lorenz, R., Takemi, T., Ito, J., and Tyler, D. (2016).  Dust Devil Formation.  Space Science Reviews, 203(1-4):183–207.
  * Rafkin and Soto, (2020) Rafkin, S. C. and Soto, A. (2020).  Air-sea interactions on titan: Lake evaporation, atmospheric circulation, and cloud formation.  Icarus, page 113903.
  * Rennó et al., (1998) Rennó, N. O., Burkett, M. L., and Larkin, M. P. (1998).  A Simple Thermodynamical Theory for Dust Devils.  Journal of Atmospheric Sciences, 55(21):3244–3252.
  * Seiff et al., (1985) Seiff, A., Schofield, J. T., Kliore, A. J., Taylor, F. W., Limaye, S. S., Revercomb, H. E., Sromovsky, L. A., Kerzhanovich, V. V., Moroz, V. I., and Marov, M. Y. (1985).  Models of the structure of the atmosphere of Venus from the surface to 100 kilometers altitude.  Advances in Space Research, 5:3–58.
  * Sergeev et al., (2020) Sergeev, D. E., Lambert, F. H., Mayne, N. J., Boutle, I. A., Manners, J., and Kohary, K. (2020).  Atmospheric convection plays a key role in the climate of tidally-locked terrestrial exoplanets: insights from high-resolution simulations.  arXiv e-prints, page arXiv:2004.03007.
  * Sinclair, (1973) Sinclair, P. C. (1973).  The Lower Structure of Dust Devils.  Journal of Atmospheric Sciences, 30(8):1599–1619.
  * Skamarock and Klemp, (2008) Skamarock, W. C. and Klemp, J. B. (2008).  A time-split nonhydrostatic atmospheric model for weather research and forecasting applications.  Journal of Computational Physics, 227:3465–3485.
  * Spiga et al., (2016) Spiga, A., Barth, E., Gu, Z., Hoffmann, F., Ito, J., Jemmett-Smith, B., Klose, M., Nishizawa, S., Raasch, S., Rafkin, S., Takemi, T., Tyler, D., and Wei, W. (2016).  Large-Eddy Simulations of Dust Devils and Convective Vortices.  Space Science Reviews, 203:245–275.
  * Spiga et al., (2010) Spiga, A., Forget, F., Lewis, S. R., and Hinson, D. P. (2010).  Structure and dynamics of the convective boundary layer on mars as inferred from large-eddy simulations and remote-sensing measurements.  Quarterly Journal of the Royal Meteorological Society, 136:414–428.
  * Spiga et al., (2011) Spiga, A., Forget, F., Madeleine, J.-B., Montabone, L., Lewis, S. R., and Millour, E. (2011).  The impact of martian mesoscale winds on surface temperature and on the determination of thermal inertia.  Icarus, 212:504–519.
  * Spiga et al., (2021) Spiga, A., Murdoch, N., Lorenz, R., Forget, F., Newman, C., Rodriguez, S., Pla-Garcia, J., Moreiras, D. V., Banfield, D., Perrin, C., Mueller, N. T., Lemmon, M., Millour, E., and Banerdt, W. B. (2021).  A Study of Daytime Convective Vortices and Turbulence in the Martian Planetary Boundary Layer Based on Half-a-Year of InSight Atmospheric Measurements and Large-Eddy Simulations.  Journal of Geophysical Research (Planets), 126(1):e06511.
  * Stull, (1988) Stull, R. B. (1988).  An introduction to boundary layer meteorology.  Springer, Dordrecht.
  * Takagi et al., (2010) Takagi, M., Suzuki, K., Sagawa, H., Baron, P., Mendrok, J., Kasai, Y., and Matsuda, Y. (2010).  Influence of CO2 line profiles on radiative and radiative-convective equilibrium states of the Venus lower atmosphere.  Journal of Geophysical Research (Planets), 115(E6):E06014.
  * Toigo et al., (2003) Toigo, A. D., Richardson, M. I., Ewald, S. P., and Gierasch, P. J. (2003).  Numerical simulation of Martian dust devils.  Journal of Geophysical Research (Planets), 108(E6):5047.
  * Weitz et al., (1994) Weitz, C. M., Plaut, J. J., Greeley, R., and Saunders, R. S. (1994).  Dunes and Microdunes on Venus: Why Were So Few Found in the Magellan Data?  Icarus, 112(1):282–295.
  * Wing et al., (2017) Wing, A. A., Emanuel, K., Holloway, C. E., and Muller, C. (2017).  Convective Self-Aggregation in Numerical Simulations: A Review.  Surveys in Geophysics, 38(6):1173–1197.
  * Wing and Emanuel, (2014) Wing, A. A. and Emanuel, K. A. (2014).  Physical mechanisms controlling self-aggregation of convection in idealized numerical modeling simulations.  Journal of Advances in Modeling Earth Systems, 6(1):59–74.
  * Yamamoto, (2011) Yamamoto, M. (2011).  Microscale simulations of Venus convective adjustment and mixing near the surface: Thermal and material transport processes.  Icarus, 211:993–1006.
  * Zhang et al., (2017) Zhang, X., Tian, F., Wang, Y., Dudhia, J., and Chen, M. (2017).  Surface Variability of Short-wavelength Radiation and Temperature on Exoplanets around M Dwarfs.  The Astrophysical Journal Letters, 837:L27.