# Distributed Data-Driven Predictive Control for Multi-Agent Collaborative
Legged Locomotion

Randall T. Fawcett1, Leila Amanzadeh1, Jeeseop Kim1, Aaron D. Ames2, and Kaveh
Akbari Hamed1 The work of R. T. Fawcett is supported by the National Science
Foundation (NSF) under Grant 2128948. The work of K. Akbari Hamed is supported
by the NSF under Grants 1924617 and 2128948. The work of A. D. Ames is
supported by the NSF under Grant 1924526.1R. T. Fawcett, L. Amanzadeh, J. Kim,
and K. Akbari Hamed (Corresponding Author) are with the Department of
Mechanical Engineering, Virginia Tech, Blacksburg, VA, 24061, USA, {randallf,
leila7, jeeseop<EMAIL_ADDRESS>D. Ames is with the Department
of Mechanical and Civil Engineering, California Institute of Technology,
Pasadena, CA, 91125, USA<EMAIL_ADDRESS>
###### Abstract

The aim of this work is to define a planner that enables robust legged
locomotion for complex multi-agent systems consisting of several holonomically
constrained quadrupeds. To this end, we employ a methodology based on
behavioral systems theory to model the sophisticated and high-dimensional
structure induced by the holonomic constraints. The resulting model is then
used in tandem with distributed control techniques such that the computational
burden is shared across agents while the coupling between agents is preserved.
Finally, this distributed model is framed in the context of a predictive
controller, resulting in a robustly stable method for trajectory planning.
This methodology is tested in simulation with up to five agents and is further
experimentally validated on three A1 quadrupedal robots subject to various
uncertainties including payloads, rough terrain, and push disturbances.

###### Index Terms:

Legged Robots, Motion Control, Multi-Contact Whole-Body Motion Planning and
Control

##  I Introduction

This work investigates multi-agent systems composed of high-dimensional
quadrupedal robots that are rigidly holonomically constrained to one another
using ball joints, introducing high interaction forces (see Fig. 1). Current
state-of-the-art approaches, even when considering only a single agent,
generally involve using a reduced-order model of some kind [1]. For a single
agent, there have been many template models that have worked effectively as
will be discussed shortly. However, template models for multi-agent systems,
particularly those with large interaction forces, have not yet been fully
explored. In particular, it is difficult to use traditional modeling
techniques to model an increasing number of agents due to increased dynamic
complexity. Even in the case that a large-scale dynamical system with strong
interaction forces could be modeled at the reduced-order level, it would
likely be of such complexity that it would no longer efficiently function as a
reduced-order model for control. The goal of this work, therefore, is to
address these issues. Namely, we aim to synthesize data-driven template models
for planning of large multi-agent systems such that the resulting planner is
computationally efficient for use in real-time.

Figure 1: Snapshot showing the locomotion of three holonomically constrained
A1 robots on wooden blocks.

###  I-A Multi-Agent Systems

Multi-Agent systems have been a major area of research, particularly in the
context of cooperative manipulation [2, 3], control of unmanned aerial
vehicles [4, 5], and ground vehicles [6, 7]. However, when considering many
agents, the system is usually high-dimensional, which also motivates
distributed approaches [8, 9, 10]. These methods reduce the computational
burden for controllers and planners, but they have also not been readily
extended to legged robotic systems with high dimensionality, unilateral
constraints, and high interaction forces, as is the case in this work.
Planning for multi-agent systems also usually necessitates the use of a
reduced-order model [11], but in the case of holonomically constrained
quadrupeds, the best model is unclear.

Figure 2: Overview of the proposed control algorithm with distributed DDPC
algorithms at the high level for trajectory optimization of cooperative
locomotion and nonlinear controllers at the low level for tracking and whole-
body motion control.

###  I-B Reduced-Order Models

In order to create predictive controllers for complex systems, it is often
necessary to utilize a reduced-order model. In terms of legged locomotion, a
very well-known and studied model is that of the linear inverted pendulum
(LIP) [12], which has been used for both quadrupedal and bipedal locomotion.
However, the LIP model generally enforces quasi-static motions due to the
requirement that the center of pressure remains within the convex hull formed
by the contacting points with the environment. Nonetheless, the model’s
simplicity has been of great interest over the years and has been used
extensively in both simulation and experiments using various platforms [13,
14, 15, 12, 16]. It is also worth mentioning that this model is not amenable
to adding moments induced about the center of mass (COM) by external forces
[17, 18], which is important when dealing with holonomically constrained
multi-agent systems.

More recently, the single rigid body (SRB) model has gained traction to
achieve dynamic locomotion [19]. This reduced-order model has proven to be
effective in creating natural and dynamic motions for various legged robot
platforms [20, 21, 22, 19]. However, a drawback is that the legs are assumed
to be massless for this model, which can become difficult when using robots
with more mass. Conversely, one major benefit of the SRB model, in contrast to
the LIP model, is that it allows one to consider the moments induced about the
COM caused by external forces, which is useful in this case because of the
holonomic constraints. Unfortunately, using the SRB model subject to
interaction forces among holonomically constrained agents becomes
prohibitively complex as the number of agents increases, which motivates the
use of data-driven methods.

###  I-C Data-Driven Methods

Data-driven methods are becoming more popular as systems become more complex
[23]. This motivates their use in multi-agent systems, especially in the case
where the individual agents are complex. Even though individual systems can be
modeled adequately, their combined dynamics make the system increasingly
sophisticated, and an effective reduced-order model has yet to be determined.
One such method recently popularized in the robotics community is
reinforcement learning [24, 25]. This method has been used for collaborative
locomotion among multiple quadrupeds [26] but takes considerable time and
computation power. Furthermore, reinforcement learning provides little
intuition as to the inner workings of the result. In this work, we pursue a
different avenue.

Here, we consider the use of data-driven methods in the context of behavioral
systems theory. This methodology parameterizes a linear time-invariant (LTI)
system directly in terms of its measured trajectories [27, 28, 29]. Behavioral
systems theory, when used in a predictive control or planning framework, is
usually referred to as data-driven predictive control (DDPC) or data-enabled
predictive control (DeePC) and has lately become of great interest to the
robotics community [30, 31, 32, 33, 34, 35, 36]. Recently, extensions have
been made to use such methods for certain types of nonlinear systems [37] and
linear parameter varying systems [38], but there have yet to be any proofs
extending to a broad class of general nonlinear systems. There have also been
advances in using these methods for stochastic implementations, and in
particular, there have been experimental validations for nonlinear systems,
even though the theory does not directly apply [31, 32, 33, 36]. However, to
the best of the authors’ knowledge, there have not been implementations for
multi-agent systems. We further note that the extension to collaborative
legged locomotion introduces many difficulties, including the hybrid nature of
legged locomotion and unilateral constraints.

###  I-D Goals, Objectives, and Contributions

The overarching goal of this paper is to develop a computationally tractable
real-time predictive planner for collaborative legged locomotion using a
behavioral approach. Namely, this work’s objectives and key contributions are
enumerated as follows: 1) A model is created using concepts from behavioral
systems theory for systems of quadrupeds that are holonomically constrained to
one another. 2) The model is used in the context of a distributed predictive
control framework such that each agent can effectively plan for its own
motions while considering the motions of other agents (see Fig. 2). 3)
Simulation results for $5$ constrained quadrupeds in the presence of
uncertainty are provided. 4) We present extensive experimental validation on a
team of $3$ holonomically constrained quadrupeds. The experimental validation
shows robust locomotion of the A1 robot subject to various uncertainties,
including rough terrain, push disturbances, and outdoor environments.

##  II Preliminaries

In this section, we overview some important concepts from behavioral systems
theory that will be used throughout this work. Behavioral systems theory
provides a manner in which data collected from a system can be leveraged to
directly create a model. In particular, consider an LTI model with the state
vector $x_{k}\in\mathbb{R}^{n}$, the input vector $u_{k}\in\mathbb{R}^{m}$,
and the output vector $y_{k}\in\mathbb{R}^{p}$ for $k\in\mathbb{Z}_{\geq
0}:=\\{0,1,\cdots\\}$. Such a model can be represented in discrete time as
follows:

$\displaystyle x_{k+1}$ $\displaystyle=A\,x_{k}+B\,u_{k}$ $\displaystyle
y_{k}$ $\displaystyle=C\,x_{k}+D\,u_{k},$ (1)

where $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times m}$,
$C\in\mathbb{R}^{p\times n}$, and $D\in\mathbb{R}^{p\times m}$ are the state-
space matrices, which are unknown. In this notation, $n$, $m$, and $p$
represent the number of states, inputs, and outputs, respectively. One
difference between traditional system identification and behavioral systems
theory is that, traditionally, one would attempt to reconstruct the matrices
of $\eqref{LTI System}$ using data. In the behavioral context, we obtain an
input-output (I-O) model without reconstructing the state matrices.

To introduce the concepts, consider $L,T\in\mathbb{Z}_{\geq 0}$, where $T\geq
L$. In addition, define some input trajectory
$u^{\textrm{d}}\in\mathbb{R}^{mT}$ composed of a sequence of
$u^{\textrm{d}}_{k}$, i.e.,
$u^{\textrm{d}}:=\textrm{col}(u^{\textrm{d}}_{0},\ldots,u^{\textrm{d}}_{T-1})$.
In this notation, “col” represents the column operator. Using this trajectory,
one can construct the following Hankel matrix

$\mathscr{H}_{L}(u^{\textrm{d}}):=\begin{bmatrix}u^{\textrm{d}}_{0}&u^{\textrm{d}}_{1}&\cdots&u^{\textrm{d}}_{T-L}\\\
u^{\textrm{d}}_{1}&u^{\textrm{d}}_{2}&\cdots&u^{\textrm{d}}_{T-L+1}\\\
\vdots&\vdots&\ddots&\vdots\\\
u^{\textrm{d}}_{L-1}&u^{\textrm{d}}_{L}&\cdots&u^{\textrm{d}}_{T-1}\end{bmatrix}\in\mathbb{R}^{mL\times{(T-L+1)}}.$
(2)

######  Definition 1 ([30])

The signal $u^{\textrm{d}}$ is said to be persistently exciting of order $L$
if $\mathscr{H}_{L}(u^{\textrm{d}})$ is full row rank.

######  Definition 2 ([30])

The sequence $\\{(u_{k},y_{k})\\}_{k=0}^{T-1}$ is said to be a trajectory of
the LTI system (II) if there exists an initial condition $x_{0}$ and a state
sequence $\\{x_{k}\\}_{k=0}^{T}$ that meets the state and output equations in
(II).

Using Definitions 1 and 2, we can now present a foundational theorem in
behavioral systems theory, used to represent an LTI system based on observed
trajectories.

###### Theorem 1

[27, Theorem 1] Let an observed trajectory of (II), referred to as data, be
denoted by $\\{(u^{\textrm{d}}_{k},y^{\textrm{d}}_{k})\\}_{k=0}^{T-1}$. If
$u^{\textrm{d}}$ is persistently exciting of order $L+n$, then
$\\{(\bar{u}_{k},\bar{y}_{k})\\}_{k=0}^{L-1}$ is a trajectory of the system if
and only if there exists $g\in\mathbb{R}^{T-L+1}$ such that

$\begin{bmatrix}\mathscr{H}_{L}(u^{\textrm{d}})\\\
\mathscr{H}_{L}(y^{\textrm{d}})\end{bmatrix}g=\begin{bmatrix}\bar{u}\\\
\bar{y}\end{bmatrix}.$ (3)

Theorem 1 provides a constructive manner in which an LTI system can be
represented through its trajectories without system identification. This
concept will be used throughout this work as we aim to parameterize a complex
system of robots by using only their trajectories. To do so, we consider two
different horizons denoted by $T_{\textrm{ini}}$ and $N$, which represent the
estimation horizon and control horizon, respectively. The estimation horizon
encapsulates the input-output pairs that are required in order to determine
the initial conditions of a trajectory given a particular I-O sequence
$\\{(\bar{u}_{k},\bar{y}_{k})\\}_{k=0}^{L-1}$ from (3). In contrast, the
prediction horizon represents how far into the future predictions of the
trajectories are made, similar to that of traditional model predictive control
(MPC). Finally, we define $L=T_{\textrm{ini}}+N$ for compact notation. We
denote the collected I-O data by $(u^{\textrm{d}},y^{\textrm{d}})$ and can
decompose the Hankel matrices of (3) into two parts as follows:

$\displaystyle\mathscr{H}_{L}(u^{\textrm{d}})=\begin{bmatrix}U_{p}\\\
U_{f}\end{bmatrix},\quad$
$\displaystyle\mathscr{H}_{L}(y^{\textrm{d}})=\begin{bmatrix}Y_{p}\\\
Y_{f}\end{bmatrix},$ (4)

where $U_{p}\in\mathbb{R}^{mT_{\textrm{ini}}\times(T-L+1)}$ and
$Y_{p}\in\mathbb{R}^{pT_{\textrm{ini}}\times(T-L+1)}$ are the portions of the
Hankel matrices used for estimating the initial condition (i.e., past), and
$U_{f}\in\mathbb{R}^{mN\times(T-L+1)}$ and
$Y_{f}\in\mathbb{R}^{pN\times(T-L+1)}$ are the portions used for prediction
(i.e., future). A necessary and sufficient condition to establish that the
Hankel matrix is sufficiently rich is to choose $T$ such that
$T\geq(m+1)(T_{\textrm{ini}}+N+n)-1$, which is a well-known result in
behavioral systems theory.

##  III Distributed DDPC for Trajectory Planning

In this section, we present our main contribution—the development of a DDPC
for constrained multi-agent systems.

###  III-A Data-Driven Predictive Control

The aim of this section is to outline how the data-driven approach is used for
predictive control, as well as highlight some difficulties. To begin, consider
the real-time DeePC methodology provided in [31, 32] as follows:

$\displaystyle\min_{(u,y,g,\sigma)}\,\,$
$\displaystyle\sum^{N-1}_{k=0}\Big{(}\|y_{k}-y^{\textrm{des}}_{k}\|^{2}_{Q}+\|u_{k}\|^{2}_{R}\Big{)}+\lambda_{g}\|g\|^{2}+\lambda_{\sigma}\|\sigma\|^{2}$
s.t. $\displaystyle\begin{bmatrix}U_{p}\\\ Y_{p}\\\ U_{f}\\\
Y_{f}\end{bmatrix}g+\begin{bmatrix}0\\\ \sigma\\\ 0\\\
0\end{bmatrix}=\begin{bmatrix}u_{\textrm{ini}}\\\ y_{\textrm{ini}}\\\ u\\\
y\end{bmatrix}$ $\displaystyle\,u_{k}\in\pazocal{U},\quad
y_{k}\in\pazocal{Y},\quad k=0,\ldots,N-1,$ (5)

where $Q\in\mathbb{R}^{p\times p}$ and $R\in\mathbb{R}^{m\times m}$ are
positive definite weighting matrices, $\|y\|^{2}_{Q}:=y^{\top}Q\,y$,
$\\{y^{\textrm{des}}_{k}\\}_{k=0}^{N-1}$ represents a desired trajectory, and
$\pazocal{U}$ and $\pazocal{Y}$ are the feasible input and output sets,
respectively. In our notation, $(u_{\textrm{ini}},y_{\textrm{ini}})$ denotes
the past measured trajectory over the estimation horizon $T_{\textrm{ini}}$,
which provides feedback directly into the model. In addition, $(u,y)$
represents the predicted I-O trajectory over the control horizon $N$. This
method has proven to be robust and has worked for several nonlinear systems
[31, 32]. One of the primary reasons for this is the addition of $\lambda_{g}$
and $\lambda_{\sigma}$, which are positive weighting factors meant to
regularize the $g$ vector from Theorem 1 and penalize the defect variable
$\sigma$, respectively. Note that $\sigma$ is added to lessen the effect that
noisy data has on the system. If the data were to contain no noise, this
variable could be removed, though it is generally required in practice.

Although this methodology has the benefit of not requiring a direct model, it
also comes with considerable computational complexity as the system increases
in size due primarily to the vector $g$. Consequently, this method is
intractable for real-time computation on teams of quadrupedal robots. For a
more in-depth discussion on the matter, we refer the interested reader to
[36]. To circumvent the problem, we adopt an offline approximation for $g$,
analogous to [32, 36] as follows:

$\displaystyle g=\begin{bmatrix}U_{p}\\\ Y_{p}\\\
U_{f}\end{bmatrix}^{\dagger}\begin{bmatrix}u_{\textrm{ini}}\\\
y_{\textrm{ini}}\\\
u\end{bmatrix},\,\,y=\pazocal{G}\begin{bmatrix}u_{\textrm{ini}}\\\
y_{\textrm{ini}}\\\ u\end{bmatrix},\,\,\pazocal{G}:$
$\displaystyle=Y_{f}\begin{bmatrix}U_{p}\\\ Y_{p}\\\
U_{f}\end{bmatrix}^{\dagger},$ (6)

where $(\cdot)^{\dagger}$ represents the pseudo inverse and $\pazocal{G}$
denotes the data-driven state transition matrix over $N$-steps. Using this
procedure, we can remove $g$ from the optimization problem (III-A),
considerably reducing the number of decision variables.

###  III-B Distributed Multi-Agent Trajectory Planning

Here, we outline how the data-driven model (6) can be used to create a
distributed planner for groups of holonomically constrained legged robots.
First, we present the control law in a centralized manner, i.e., assuming that
one planner is used to control the whole system. This is later decomposed for
distributed computation. In particular, consider the centralized predictive
control problem utilizing (6) as follows:

$\displaystyle\min_{(u,y)}\quad$
$\displaystyle\sum^{N-1}_{k=0}\Big{(}\|y_{k}-y_{k}^{\textrm{des}}\|^{2}_{Q}+\|u_{k}\|^{2}_{R}\Big{)}$
s.t. $\displaystyle y=\pazocal{G}\begin{bmatrix}u_{\textrm{ini}}\\\
y_{\textrm{ini}}\\\ u\end{bmatrix}$
$\displaystyle\,u_{k}\in\pazocal{U},\,\,y_{k}\in\pazocal{Y},\quad\,\,k=0,\ldots,N-1.$
(7)

This method has been shown to be amenable to trajectory planning for single-
agent legged robots [36] but has not been used for multi-agent systems. In the
multi-agent context, the state transition matrix $\pazocal{G}$ is created
using data from all of the agents. In particular, Theorem 1 considers the
trajectory $(u^{\textrm{d}},y^{\textrm{d}})$ to define the Hankel matrix. To
construct the Hankel martix for multi-agent systems, we can define
$u^{\textrm{d}}:=\\{\textrm{col}(u^{\textrm{d},1}_{k},u^{\textrm{d},2}_{k},\cdots,u^{\textrm{d},n_{a}}_{k})\\}_{k=0}^{T-1}$
and
$y^{\textrm{d}}:=\\{\textrm{col}(y^{\textrm{d},1}_{k},y^{\textrm{d},2}_{k},\cdots,y^{\textrm{d},n_{a}}_{k})\\}_{k=0}^{T-1}$,
where $(\cdot)^{\textrm{d},i}_{k}$ for all
$i\in\pazocal{I}:=\\{1,\cdots,n_{a}\\}$ denotes the data contributed by agent
$i$ at the sample time $k$, and $n_{a}$ is the total number of agents. Using
these combined I-O pairs, we obtain a large Hankel matrix describing the
entire complex system. Furthermore, this new Hankel matrix can be decomposed
according to (4), and the corresponding $g$ vector can be approximated using
(6).

In moving to multi-agent systems, however, it is desirable to share the
computational load between agents. In order to do so, we consider a
decomposition of $\pazocal{G}$ as follows:

$\displaystyle\pazocal{G}=\begin{bmatrix}\pazocal{G}_{1,1}&\pazocal{G}_{1,2}&\cdots&\pazocal{G}_{1,n_{a}}\\\
\pazocal{G}_{2,1}&\pazocal{G}_{2,2}&\cdots&\pazocal{G}_{2,n_{a}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\pazocal{G}_{n_{a},1}&\pazocal{G}_{n_{a},2}&\cdots&\pazocal{G}_{n_{a},n_{a}}\end{bmatrix},$
(8)

where $\pazocal{G}_{i,j}$ represents the effect of agent $j$ on the predicted
output of agent $i$ (i.e., coupling). In particular, the predicted output of
agent $i\in\pazocal{I}$ can be expressed by

$y^{i}=\pazocal{G}_{i,i}\begin{bmatrix}u^{i}_{\textrm{ini}}\\\
y^{i}_{\textrm{ini}}\\\ u^{i}\end{bmatrix}+\sum_{j\neq
i}\pazocal{G}_{i,j}\begin{bmatrix}u^{j}_{\textrm{ini}}\\\
y^{j}_{\textrm{ini}}\\\ u^{j}\end{bmatrix},$ (9)

where $(\cdot)^{i}$ represents the variables related to the $i$th agent. To
alleviate the coupling problem, we adopt a one-step communication delay
protocol. In particular, we assume that at every time sample $t$, each local
DDPC has access to the optimal predicted solutions of the other local DDPCs
and their past measurements at time $t-1$. Using this assumption, the
predicted output of agent $i$ can be approximated by

$y^{i}\approx\pazocal{G}_{i,i}\begin{bmatrix}u^{i}_{\textrm{ini}}\\\
y^{i}_{\textrm{ini}}\\\ u^{i}\end{bmatrix}+\sum_{j\neq
i}\pazocal{G}_{i,j}\begin{bmatrix}u^{j}_{\textrm{ini}|t-1}\\\
y^{j}_{\textrm{ini}|t-1}\\\ u^{j}_{|t-1}\end{bmatrix},$ (10)

where $u^{j}_{|t-1}$ denotes the optimal solution of the local DDPC for agent
$j$ at time $t-1$, wherein the summation is constant until the next time the
DDPC is run. In addition,
$(u^{j}_{\textrm{ini}|t-1},y^{j}_{\textrm{ini}|t-1})$ represents the past
measurements of agent $j$ at time $t-1$. This choice allows each local planner
to run independently and results in the following network of distributed DDPCs
(local QPs) to be solved at time $t$

$\displaystyle\min_{(u^{i},y^{i})}\quad$
$\displaystyle\sum^{N-1}_{k=0}\Big{(}\|y^{i}_{k}-y_{k}^{\textrm{des},i}\|^{2}_{Q}+\|u^{i}_{k}-u_{k}^{\textrm{des},i}\|^{2}_{R}\Big{)}$
s.t. $\displaystyle
y^{i}=\pazocal{G}_{i,i}\begin{bmatrix}u^{i}_{\textrm{ini}}\\\
y^{i}_{\textrm{ini}}\\\ u^{i}\end{bmatrix}+\sum_{j\neq
i}\pazocal{G}_{i,j}\begin{bmatrix}u^{j}_{\textrm{ini}|t-1}\\\
y^{j}_{\textrm{ini}|t-1}\\\ u^{j}_{|t-1}\end{bmatrix}$ $\displaystyle
u^{i}_{k}\in\pazocal{U},\,\,y^{i}_{k}\in\pazocal{Y},\quad k=0,\ldots,N-1,$
(11)

where $u_{k}^{\textrm{des},i}$ and $y_{k}^{\textrm{des},i}$ represent the
desired inputs and outputs for agent $i$ at the prediction step $k$ (see Fig.
2). The optimal input and output trajectories are then passed to the low-level
controller for tracking. In this work, we choose a subset of COM state
variables for the output $y$ while taking the ground reaction forces (GRFs) as
the control inputs $u$. This will be clarified more in Section V-A.
Consequently, the feasible set $\pazocal{U}$ is chosen as the linearized
friction cone, i.e.,
$\pazocal{U}=\pazocal{FC}:=\\{\textrm{col}(f_{x},f_{y},f_{z})|f_{z}>0,\,|f_{x}|\leq\frac{\mu}{\sqrt{2}}f_{z},\,|f_{y}|\leq\frac{\mu}{\sqrt{2}}f_{z}\\}$,
where $\mu$ denotes the friction coefficient.

Figure 3: (a) Simulation results of $5$ agents over varying terrain with a
payload of $10$ (kg), (b) rough terrain experiment with unstructured wooden
blocks and a $4.5$ (kg) payload, (c) experiment with push disturbances, and
(d) experiment maneuvering over gravel. Videos are available online [39].

##  IV Nonlinear Low-Level Controller

The purpose of this section is to briefly provide the details of the low-level
controller [40], based on QP and virtual constraints [41], to be used in this
work. In particular, we model each robot using a floating base and represent
the generalized coordinates of the system by
$q\in\pazocal{Q}\subset\mathbb{R}^{n_{q}}$, where $\pazocal{Q}$ is the
configuration space and $n_{q}$ represents the number of degrees of freedom
(DOFs) of the system. We further denote the joint-level torques by
$\tau\in\pazocal{T}\subset\mathbb{R}^{m_{\tau}}$. In this notation,
$\pazocal{T}$ denotes the allowable torques and $m_{\tau}$ denotes the number
of actuators. It is assumed that the interaction wrenches between agents are
encapsulated by the data-driven model. Therefore, we neglect these wrenches in
the low-level controller. In particular, the overarching equations of motion
for the full-order system of each agent become

$M(q)\,\ddot{q}+H(q,\dot{q})=\Upsilon\,\tau+J^{\top}(q)\,f,$ (12)

where $M(q)\in\mathbb{R}^{n_{q}\times n_{q}}$ is the mass-inertia matrix,
$H(q,\dot{q})\in\mathbb{R}^{n_{q}}$ denotes the Coriolis, centrifugal, and
gravitational terms, $\Upsilon\in\mathbb{R}^{n_{q}\times{}m_{\tau}}$
represents the input matrix, $J(q)$ denotes the contact Jacobian matrix, and
$f$ represents the GRFs at the stance feet. We further suppose that the
positions of the stance leg ends, denoted by $r$, do not slip, i.e.,
$\ddot{r}=0$.

With the dynamics in hand, we can now present the virtual constraints
controller. In particular, we aim to track both force and COM trajectories
generated by the high-level distributed planners. For this purpose, we
consider virtual constraints as output functions to be regulated as
$h(q,t):=h_{0}(q)-h^{\textrm{des}}(t)$. These virtual constraints are then
imposed via partial feedback linearization [42]. Here, $h_{0}(q)$ denotes the
controlled variables consisting of the COM position, orientation, and the
Cartesian coordinates of the swing leg ends. Finally, $h^{\textrm{des}}(t)$
represents the desired evolution of $h_{0}(q)$. In this work, the desired end
position for a swing leg is chosen using the Raibert heuristic [43, Eq. (4),
pp. 46], and the trajectory for the swing leg is defined using a Bézier
polynomial. These virtual constraints, along with the no slippage condition,
are concatenated into a single strictly convex QP to be solved at 1kHz as
follows [40]

$\displaystyle\min_{(\tau,f,\delta)}\,\,\,$
$\displaystyle\frac{\gamma_{1}}{2}\|\tau\|^{2}+\frac{\gamma_{2}}{2}\|f-f^{\textrm{des}}\|^{2}+\frac{\gamma_{3}}{2}$
$\displaystyle\|\delta\|^{2}$ s.t.
$\displaystyle\ddot{h}(\tau,f)=-K_{P}\,h-K_{D}\,\dot{h}+\delta$ (Output
Dynamics) $\displaystyle\ddot{r}(\tau,f)=0$ (No slippage)
$\displaystyle\tau\in\pazocal{T},\quad f\in\pazocal{FC}$
$\displaystyle\textrm{(Feasibility)},$ (13)

where $\gamma_{1}$, $\gamma_{2}$, and $\gamma_{3}$ are positive weighting
factors, and the desired force profile $f^{\textrm{des}}(t)$ (i.e., inputs
$u$) is prescribed by the high-level DDPC in (III-B). The equality constraints
are expressed as 1) the output dynamics
$\ddot{h}+K_{D}\,\dot{h}+K_{P}\,h=\delta$ for positive definite gain matrices
$K_{P}$ and $K_{D}$ and $\delta$ being a defect variable to ensure
feasibility, and 2) the no slippage condition $\ddot{r}=0$. We remark that
$\ddot{h}$ and $\ddot{r}$ are affine functions of $(\tau,f)$, hence, the
problem is convex. We direct the reader to [17, Appendix A] for more
information regarding the derivation of $\ddot{h}$ and $\ddot{r}$ according to
Lie derivatives. The QP solves for the minimum-power torques $\tau$ while
tracking the prescribed forces and COM trajectory.

##  V Experiments

In this section, we provide the procedure for collecting the data for the
reduced-order model, and further, provide the simulation and experimental
results. Here we consider the $18$-DOF quadruped A1 made by Unitree. The robot
is modeled using a floating base, with the first $6$ DOFs being composed of
the unactuated position and orientation of the trunk. The remaining DOFs are
composed of the actuated hip roll, hip pitch, and knee pitch joints for each
leg. The robot weighs around $12.45$ (kg) and the center of the trunk is
$0.26$ (m) above the ground during locomotion. We are, however, interested in
multi-agent systems. We assume that each of the agents is holonomically
constrained to other agents in a complete graph, i.e., there is no relative
translational motion between agents. In order to accomplish this, we connect
the robots together rigidly through a ball joint (see Fig. 3). In simulations,
this is imposed via a distance constraint.

###  V-A Data Collection

We begin by describing the I-O pairs considered in this work and the manner in
which the data was collected. The data that is to be used in the Hankel matrix
is collected in simulation at $100$ (Hz) during nominal locomotion. This is in
contrast to [36], which considered data collected on hardware. However, one of
the contributions of this work is to create a reduced-order model of highly
complex systems that interact with one another. From this standpoint, it makes
sense to consider simulation data since each individual agent can be modeled
accurately, but the complexity of collaborating agents is prohibitive in terms
of defining a physics-based reduced-order model. This also shows good sim-to-
real transfer for the learned model, as will be shown in what follows.

In order to perform the data collection in simulation, we first choose the
inputs to be the forces at the contacting leg ends and and we take the outputs
as
$y^{\textrm{d}}=\textrm{col}(\textrm{z},\dot{\textrm{x}},\dot{\textrm{y}},\textrm{roll},\textrm{pitch},\omega_{z})$,
where z is the standing height, $\dot{\textrm{x}}$ and $\dot{\textrm{y}}$ are
the linear velocities of the COM in the transverse plane, and $\omega_{z}$
denotes the angular velocity about the vertical axis of the torso. These
outputs are chosen because they represent the variables of interest to an end-
user when providing joystick commands to a robot. It should be noted that
other I-O realization could provide fruitful results depending on the goal of
the planner and the available measurements. During the data collection
procedure, the quadrupeds are commanded to walk around in RaiSim [44] using
just the low-level controller (IV), while random noise is injected into the
desired forces, which helps ensure the persistence of excitation. In
particular, we choose the desired forces to be
$u^{\textrm{des},i}_{k,\ell}:=\textrm{col}(0,0,\frac{m^{\textrm{net},i}g_{0}}{N^{i}_{c,k}})$
for each contacting leg $\ell\in\pazocal{C}_{k}^{i}$ and zero otherwise, where
$m^{\textrm{net},i}$ is the total mass of agent $i$, $g_{0}$ is the
gravitational acceleration, $N_{c,k}^{i}$ is the anticipated number of
contacting legs for agent $i$ at time $k$, and $\pazocal{C}_{k}^{i}$ is the
anticipated set of contacting legs for agent $i$ at time $k$. That is, when
collecting data, we choose the desired force to be a random perturbation about
the nominal amount of force that is required to hold the quadruped in a static
position based on the number of anticipated contacts.

For this problem, we consider an estimation horizon of $T_{\textrm{ini}}=10$
and $N=25$ for the prediction horizon. We further ensure that $T$ is chosen
such that the amount of data collected far exceeds that required by the
general theory. This, in turn, assists in potentially capturing more nonlinear
information while it also reduces the impact of noise. The collection of
additional data when utilizing the formulation (III-A) could pose an issue due
to an increase in decision variables. However, this is mitigated by utilizing
the approximation (6).

Figure 4: The trajectory from the planner (a) and the vertical GRF of the
front right leg from the planner (b) for agent $1$. A forward speed of $0.5$
(m/s) is commanded, and the standing height is $0.26$ (m). The multi-agent
system is subject to uneven terrain and a payload of $10$ (kg), and can
maneuver robustly. A snapshot of the simulation can be found in Fig. 3 (a).

###  V-B Simulation Experiments

The high-level planner contains $450$ decision variables per agent. In the
centralized case using $5$ agents shown here, that amounts to $2250$ decision
variables, which further motivates the necessity for a distributed approach.
Under the distributed scheme, the predictive controller is updated every $40$
(ms) ($25$ (Hz)), and the first $4$ time steps are implemented, i.e., we
predict over a horizon of $250$ (ms) and implement the first $40$ (ms) of the
prediction. The high-level planner is solved using OSQP [45] and takes
approximately $15$ (ms) on an external laptop with an Intel® Core™ i7-1185G7
running at $3.00$ GHz and $16$ GB of RAM. However, solve times of $\sim 30$
(ms) have been observed, further motivating the decision to update the planner
at $40$ (ms). The predictive controller parameters are chosen as
$Q=\textrm{diag}(1e6,1e5,1e5,2e5,1e5,1e4)$ and
$R=I\otimes\textrm{diag}(0.05,0.05,0.5)$, where $I$ is the identity of
appropriate size, and $\otimes$ represents the Kronecker product. Finally, the
parameters used by the low-level controller to track the trajectory and forces
from the planner are chosen to be $\gamma_{1}=10^{2}$, $\gamma_{2}=10^{3}$,
and $\gamma_{3}=10^{6}$, resulting in stable locomotion. All of the gains used
in simulation are the same as those used on hardware.

In order to show the efficacy of the proposed controller for $5$ agents, we
consider a compound experiment in simulation such that the multi-agent system
is subject to an unknown payload of $10$ (kg) and uneven terrain. A snapshot
of the simulation can be found in Fig. 3 (a), while the prescribed forces and
trajectory for the first agent can be found in Fig. 4. From these figures, it
is evident that the planner produces forces that are feasible, while also
resulting in a viable COM trajectory for the low-level controller to track.

Figure 5: The trajectory from the planner of agent $1$ while trotting at
approximately $0.4$ (m/s) subject to (a) rough terrain with unstructured
wooden blocks and (b) loose gravel. We further show the vertical GRF for the
front right leg produced by the planner for rough terrain in (c). The GRF
during the gravel experiment is similar. Figure 6: The plot shows the tracking
performance of the planner when using a time-varying trajectory provided using
a joystick. In this experiment, the quadrupeds navigate flat ground subject to
a $6.8$ (kg) payload. The plot shows the trajectory of agent $1$.

###  V-C Hardware Experiments

Finally, we provide hardware experiments to show the effectiveness of the
planner. In particular, for hardware experiments, we consider the use of $3$
quadrupeds that are holonomically constrained using a ball joint. Snapshots of
several experiments can be found in Fig. 3 (b)-(d), which shows the multi-
agent system subject to external disturbances and unknown environments. In
Fig. 5, we illustrate the behavior of the planner in terms of its prescribed
trajectory when walking over unstructured wooden blocks (Fig. 5 (a)) and
navigating over loose gravel (Fig. 5 (b)). The corresponding forces for the
rough terrain experiment can be found in Fig. 5 (c). Finally, we provide an
additional experiment to show the ability of the planner to track a time-
varying trajectory subject to a $6.8$ (kg) payload, where the trajectory
produced by the planner can be found in Fig. 6. It is evident that the planner
provides a robustly stable output even in the presence of significant
uncertainty in terms of payloads and various environmental factors. Videos of
the simulation and hardware experiments can be found in [39].

##  VI Conclusion

This work presented a data-driven planner for robust multi-agent quadrupedal
locomotion, wherein the robots were constrained to one another with a ball
joint. We considered the use of behavioral systems theory to model the complex
system and further proposed a distributed scheme to spread the computational
load. We provided extensive experiments both in simulation and on hardware,
which showed the robustness of this method to uncertainty in terrain,
payloads, and external disturbances. Future work will examine this methodology
when the robots do not form a complete graph, i.e., the agents are constrained
but have limited freedom to change formation. Additionally, we will explore
how this method could extend to an even greater number of agents.

## References

  * [1] R. J. Full and D. E. Koditschek, “Templates and anchors: Neuromechanical hypotheses of legged locomotion on land,” _Journal of Experimental Biology_ , vol. 202, pp. 3325–3332, December 1999.
  * [2] P. Culbertson, J.-J. Slotine, and M. Schwager, “Decentralized adaptive control for collaborative manipulation of rigid bodies,” _IEEE Transactions on Robotics_ , vol. 37, no. 6, pp. 1906–1920, 2021.
  * [3] R. Murray, Z. Li, and S. S.S., _A Mathematical Introduction to Robotic Manipulation_. Taylor & Francis/CRC, 1994\.
  * [4] Y. Chen, A. Singletary, and A. D. Ames, “Guaranteed obstacle avoidance for multi-robot operations with limited actuation: A control barrier function approach,” _IEEE Control Systems Letters_ , vol. 5, no. 1, pp. 127–132, 2020\.
  * [5] K. Sreenath and V. Kumar, “Dynamics, control and planning for cooperative manipulation of payloads suspended by cables from multiple quadrotor robots,” in _Robotics: Science and Systems (RSS)_ , 2013.
  * [6] D. Panagou, M. Turpin, and V. Kumar, “Decentralized goal assignment and trajectory generation in multi-robot networks: A multiple lyapunov functions approach,” in _IEEE International Conference on Robotics and Automation_ , Hong Kong, China, June 2014, pp. 6757–6762.
  * [7] T. Machado, T. Malheiro, S. Monteiro, W. Erlhagen, and E. Bicho, “Multi-constrained joint transportation tasks by teams of autonomous mobile robots using a dynamical systems approach,” in _IEEE International Conference on Robotics and Automation_ , 2016, pp. 3111–3117.
  * [8] W. B. Dunbar and R. M. Murray, “Distributed receding horizon control for multi-vehicle formation stabilization,” _Automatica_ , vol. 42, no. 4, pp. 549–558, 2006.
  * [9] F. Bullo, J. Cortés, and M. S, _Distributed Control of Robotic Networks: A Mathematical Approach to Motion Coordination Algorithms_. Princeton University Press, 2009.
  * [10] M. Mesbahi and E. M, _Graph Theoretic Methods in Multiagent Networks_. Princeton University Press, 2010.
  * [11] A. Perizzato, M. Farina, and R. Scattolini, “Formation control and collision avoidance of unicycle robots with distributed predictive control,” _IFAC-PapersOnLine_ , vol. 48, no. 23, pp. 260–265, 2015.
  * [12] S. Kajita, F. Kanehiro, K. Kaneko, K. Fujiwara, K. Harada, K. Yokoi, and H. Hirukawa, “Biped walking pattern generation by using preview control of zero-moment point,” in _IEEE International Conference on Robotics and Automation_ , vol. 2, Sep. 2003, pp. 1620–1626.
  * [13] R. J. Griffin, G. Wiedebach, S. Bertrand, A. Leonessa, and J. Pratt, “Walking stabilization using step timing and location adjustment on the humanoid robot, Atlas,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , Sep. 2017, pp. 667–673.
  * [14] J. Pratt, J. Carff, S. Drakunov, and A. Goswami, “Capture point: A step toward humanoid push recovery,” in _IEEE-RAS International Conference on Humanoid Robots_ , Dec 2006, pp. 200–207.
  * [15] J. Englsberger, C. Ott, M. A. Roa, A. Albu-Schäffer, and G. Hirzinger, “Bipedal walking control based on capture point dynamics,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , Sep. 2011, pp. 4420–4427.
  * [16] K. Akbari Hamed, J. Kim, and A. Pandala, “Quadrupedal locomotion via event-based predictive control and QP-based virtual constraints,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 3, pp. 4463–4470, 2020\.
  * [17] J. Kim and K. Akbari Hamed, “Cooperative locomotion via supervisory predictive control and distributed nonlinear controllers,” _Journal of Dynamic Systems, Measurement, and Control_ , vol. 144, no. 3, 2021.
  * [18] R. T. Fawcett, A. Pandala, J. Kim, and K. Akbari Hamed, “Real-Time Planning and Nonlinear Control for Quadrupedal Locomotion With Articulated Tails,” _Journal of Dynamic Systems, Measurement, and Control_ , vol. 143, no. 7, February 2021.
  * [19] G. Bledt, M. J. Powell, B. Katz, J. Di Carlo, P. M. Wensing, and S. Kim, “MIT Cheetah 3: Design and control of a robust, dynamic quadruped robot,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , Oct 2018, pp. 2245–2252.
  * [20] R. Grandia, A. J. Taylor, A. D. Ames, and M. Hutter, “Multi-layered safety for legged robots via control barrier functions and model predictive control,” in _IEEE International Conference on Robotics and Automation_ , 2021, pp. 8352–8358.
  * [21] M. Chignoli and P. M. Wensing, “Variational-based optimal control of underactuated balancing for dynamic quadrupeds,” _IEEE Access_ , vol. 8, pp. 49 785–49 797, 2020.
  * [22] Y. Ding, A. Pandala, C. Li, Y.-H. Shin, and H.-W. Park, “Representation-free model predictive control for dynamic motions in quadrupeds,” _IEEE Transactions on Robotics_ , vol. 37, no. 4, pp. 1154–1171, 2021.
  * [23] Z.-S. Hou and Z. Wang, “From model-based control to data-driven control: Survey, classification and perspective,” _Information Sciences_ , vol. 235, pp. 3–35, 2013.
  * [24] J. Lee, J. Hwangbo, L. Wellhausen, V. Koltun, and M. Hutter, “Learning quadrupedal locomotion over challenging terrain,” _Science Robotics_ , vol. 5, no. 47, 2020.
  * [25] T. Miki, J. Lee, J. Hwangbo, L. Wellhausen, V. Koltun, and M. Hutter, “Learning robust perceptive locomotion for quadrupedal robots in the wild,” _Science Robotics_ , vol. 7, no. 62, 2022.
  * [26] Y. Ji, B. Zhang, and K. Sreenath, “Reinforcement learning for collaborative quadrupedal manipulation of a payload over challenging terrain,” in _IEEE International Conference on Automation Science and Engineering_ , 2021, pp. 899–904.
  * [27] J. C. Willems, P. Rapisarda, I. Markovsky, and B. L. De Moor, “A note on persistency of excitation,” _Systems & Control Letters_, vol. 54, no. 4, pp. 325–329, 2005.
  * [28] J. C. Willems, “From time series to linear system—Part I. Finite dimensional linear time invariant systems,” _Automatica_ , vol. 22, no. 5, pp. 561–580, 1986.
  * [29] I. Markovsky, J. C. Willems, S. Van Huffel, and B. De Moor, _Exact and approximate modeling of linear systems: A behavioral approach_. Society for Industrial and Applied Mathematics, 2006\.
  * [30] J. Berberich, J. Köhler, M. A. Müller, and F. Allgöwer, “Data-driven model predictive control with stability and robustness guarantees,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 4, pp. 1702–1717, 2020.
  * [31] J. Coulson, J. Lygeros, and F. Dörfler, “Data-enabled predictive control: In the shallows of the DeePC,” in _European Control Conference_ , 2019, pp. 307–312.
  * [32] L. Huang, J. Coulson, J. Lygeros, and F. Dörfler, “Data-enabled predictive control for grid-connected power converters,” in _IEEE Conference on Decision and Control_ , 2019, pp. 8130–8135.
  * [33] J. Coulson, J. Lygeros, and F. Dörfler, “Distributionally robust chance constrained data-enabled predictive control,” _IEEE Transactions on Automatic Control_ , 2021.
  * [34] H. J. Van Waarde, J. Eising, H. L. Trentelman, and M. K. Camlibel, “Data informativity: A new perspective on data-driven analysis and control,” _IEEE Transactions on Automatic Control_ , vol. 65, no. 11, pp. 4753–4768, 2020.
  * [35] L. Wei, Y. Yan, and J. Bao, “A data-driven predictive control structure in the behavioral framework,” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 152–157, 2020.
  * [36] R. T. Fawcett, K. Afsari, A. D. Ames, and K. A. Hamed, “Toward a data-driven template model for quadrupedal locomotion,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 3, pp. 7636–7643, 2022.
  * [37] J. Berberich and F. Allgöwer, “A trajectory-based framework for data-driven system analysis and control,” in _European Control Conference_ , 2020, pp. 1365–1370.
  * [38] C. Verhoek, H. S. Abbas, R. Tóth, and S. Haesaert, “Data-driven predictive control for linear parameter-varying systems,” _IFAC-PapersOnLine_ , vol. 54, no. 8, pp. 101–108, 2021.
  * [39] Distributed data-driven predictive control for multi-agent collaborative legged locomotion. [Online]. Available: https://youtu.be/k1VsgvbR1Rs.
  * [40] R. T. Fawcett, A. Pandala, A. D. Ames, and K. Akbari Hamed, “Robust stabilization of periodic gaits for quadrupedal locomotion via QP-based virtual constraint controllers,” _IEEE Control Systems Letters_ , vol. 6, pp. 1736–1741, 2021.
  * [41] E. Westervelt, J. Grizzle, C. Chevallereau, J. Choi, and B. Morris, _Feedback Control of Dynamic Bipedal Robot Locomotion_. Taylor & Francis/CRC, 2007.
  * [42] A. Isidori, _Nonlinear Control Systems_. Springer; 3rd edition, 1995.
  * [43] M. H. Raibert, _Legged robots that balance_. MIT press, 1986.
  * [44] J. Hwangbo, J. Lee, and M. Hutter, “Per-contact iteration method for solving contact dynamics,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 2, pp. 895–902, April 2018.
  * [45] B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd, “OSQP: an operator splitting solver for quadratic programs,” _Mathematical Programming Computation_ , vol. 12, no. 4, pp. 637–672, 2020.


# Gas flow in Martian spider formation

Nicholas Attree111Corresponding author<EMAIL_ADDRESS>Erika Kaufmann
Axel Hagermann  Faculty of Natural Sciences, University of Stirling, UK  Now
at Institute for Space Research Graz, Austrian Academy of Sciences, Austria
Luleå University of Technology, Space Campus Kiruna, Sweden

###### Abstract

Martian araneiform terrain, located in the Southern polar regions, consists of
features with central pits and radial troughs which are thought to be
associated with the solid state greenhouse effect under a CO2 ice sheet.
Sublimation at the base of this ice leads to gas buildup, fracturing of the
ice and the flow of gas and entrained regolith out of vents and onto the
surface. There are two possible pathways for the gas: through the gap between
the ice slab and the underlying regolith, as proposed by Kieffer [2007], or
through the pores of a permeable regolith layer, which would imply that
regolith properties can control the spacing between adjacent spiders, as
suggested by Hao et al. [2019]. We test this hypothesis quantitatively in
order to place constraints on the regolith properties. Based on previously
estimated flow rates and thermophysical arguments, we suggest that there is
insufficient depth of porous regolith to support the full gas flow through the
regolith. By contrast, free gas flow through a regolith–ice gap is capable of
supplying the likely flow rates for gap sizes on the order of a centimetre.
This size of gap can be opened in the centre of a spider feature by gas
pressure bending the overlying ice slab upwards, or by levitating it entirely
as suggested in the original Kieffer [2007] model. Our calculations therefore
support at least some of the gas flowing through a gap opened between the
regolith and ice. Regolith properties most likely still play a role in the
evolution of spider morphology, by regolith cohesion controlling the erosion
of the central pit and troughs, for example.

###### keywords:

Mars; Mars, surface; Mars, polar geology: ices

††journal: Icarus

##  1 Introduction

Araneiform terrain consists of groups of so–called spiders; radial, dendritic
arrangements of troughs connecting to central pits, found at southern high
latitudes of Mars [Kieffer, 2000]. Their occurence was thought to be limited
to the South Polar Layered Deposits (SPLD). Spiders were also discovered
outside the area of the SPLD, although the extent of the SPLD was subsequently
re–defined. Tanaka et al. [2014] or Schwamb et al. [2018] provide more insight
the geography of spiders and the SPLD. These features are perennial, visible
all year round, scoured up to several metres deep into the regolith, and can
be tens to hundreds of metres in size [Hansen et al., 2010]. A slab of CO2 ice
begins to form in the autumn and covers this terrain during the southern
winter and spring [Kieffer, 2000]. Based on this hypothesis, Piqueux et al.
[2003] were the first to define and map the distribution of the spider
features, confirming the correlation of spiders with the existence of highly
transparent CO2 slab ice, which overlies the largely unconsolidated
particulate regolith of the South Polar Layered Deposits.

Over the years, what is now commonly referred to as the Kieffer model has
undergone extensive extensions and refinements, but there is general consensus
that spiders are formed by gas flowing under the CO2 ice [Piqueux et al.,
2003, Kieffer et al., 2006, Kieffer, 2007, Piqueux and Christensen, 2008,
Portyankina et al., 2010, Thomas et al., 2010, Pilorget et al., 2011, Thomas
et al., 2011b, Martínez et al., 2012, Pilorget and Forget, 2016]. Sub–ice gas
production is enabled by the transparency of $CO_{2}$ ice (see e.g. Hansen
[2005] and references therein). The ice slab allows spring–time insolation to
penetrate to its base and warm it in a solid-state greenhouse, a phenomenon
similar to the atmospheric greenhouse effect [Brown and Matson, 1987, Fanale
et al., 1990]. This means that sunlight can penetrate the CO2 ice layer down
to a dust deposit at depth within or below the ice, where the radiation is
absorbed. The temperature increase here leads to sublimation of the CO2 ice on
the boundary between the ice layer and the underlying material. Sublimating
gas then builds up until the pressure fractures the ice, opening a vent to the
surface, followed by gas flow from the surrounding area into the vent and
ejection. Regolith particles become entrained in the flow, being ejected along
with the gas, whilst also eroding the troughs and pits [Thomas et al., 2010].
Larger and heavier grains end up forming dark fans on top of the ice. These
can be more radial or unidirectional, depending on local surface winds
carrying the material over the top of the ice sheet. whilst the smallest and
lightest grains are ejected into the atmosphere [Portyankina et al., 2010].
This process is thought to repeat seasonally, so spiders play an important
role in Mars’ global dust cycle and therefore the planet’s climate (see also
[Kieffer, 2000, Piqueux et al., 2003]). Over time, spiders can effect
substantial morphological change and they are representative of a type
$CO_{2}$–sublimation associated erosion representative for Mars. In fact, most
morphological changes on Mars are considered to be $CO_{2}$ related today
[Piqueux and Christensen, 2008, Hansen, 2005].

As for the gas trapped unerneath the ice slab, there are several escape routes
to consider. In the commonly accepted Kieffer model 2007, gas pressure
physically lifts the ice slab, levitating at least part of the overlying slab
and allowing gas to flow through an ice–regolith gap. This allows a radially
converging flow to a central vent to drain a large area of sublimating ice
through only a small ($\sim$cm thick) gap. However, we need to bear in mind
the porous nature of Mars’ regolith, which might permit gas to flow within the
substrate underlying the ice. The importance of the latter gas escape route
has recently been stressed by Hao et al. [2019], who established a link
between spider spatial density and regolith morphology. They observed that
spiders are grouped at mutual distances smaller than that of random spacing,
suggesting a controlling length–scale whereby a spider inhibits the formation
of new spiders out to some distance from itself. Since this distance was seen
to vary between different classes of spider morphology, Hao et al. [2019]
proposed that it is controlled by the regolith properties, e.g. that different
regolith porosities/permeabilities constrain the gas flow by different
amounts, and therefore the radius around each spider where gas can build up.
Evidence for the relevance of regolith morphology was observed by Hao et al.
[2020]. These recent results highlight that the potential effectiveness of gas
flow through the regolith itself — rather than through a regolith–ice gap —
has to be addressed because it seems vital for our understanding of spiders.

We approach the problem by providing suitable limits on the effectiveness of
the underground flow of gas, based on constraints on regolith properties such
as porosity/permeability combined with some known gas physics. Regolith
porosity is a function of depth and therefore our treatment of the problem has
to involve some stratigraphic considerations. We derive the depth of porous
regolith needed to support the magnitudes of flow rate previously calculated
in hydrodynamics simulations [Thomas et al., 2011a, b] and energy balance
models [Kieffer, 2007, Pilorget et al., 2011, Pilorget and Forget, 2016]. We
then compare this to the required gap size to support laminar free–flow
between the ice and regolith, and to the expected bending of the ice sheet in
response to the rise in gas pressure beneath. We present the calculations in
Section 2, before discussing the results and their implications in Section 3,
and concluding in Section 4.

##  2 Modelling

The most common model for spider formation is illustrated in Figure 1; CO2
slab ice condenses on-top of regolith over the Martian winter. When insolation
begins to penetrate this transparent ice layer in spring, CO2 ice at the base
of the slab and top of the regolith starts to sublimate, raising the gas
pressure to the point where it exceeds the strength of the ice slab and
fracturing occurs. In the model of Kieffer [2007], this gas pressure levitates
the slab above the regolith surface by a distance of $d_{gap}$, estimated at
around one or two centimetres. Gas can also diffuse into the permeable part of
the regolith [Hao et al., 2019], which extends down to some depth $d$ (also on
the order of centimetres [Kieffer, 2007, Portyankina et al., 2010]), below
which the regolith is made impermeable by water ice filling the pore spaces.
See Mellon et al. [2004] for a detailed discussion of the geography of ground
ice on Mars. Irrespective of whether the gas flows above or through the
regolith to begin with, once fracturing has opened a vent to the surface, a
pressure gradient across the radius of the spider drives gas flow to it,
either as free–flow through the gap or viscous flow through the porous
regolith. In the latter case, flow is limited by the pressure gradient and the
regolith permeability, and we estimate these below.

Figure 1: Conceptual model for spider formation based on a layer of
transparent CO2 ice of thickness $z$

, overlaying permeable regolith of thickness $d$, above an impermeable
ice–rich regolith layer shaded dark. Gas pressure in the regolith increases
from near–atmospheric pressure around the spider vent to a maximum value at
the edge of the spider. Some models, such as those by Hao et al. [2019] and
Pilorget et al. [2011] assume no gap between ice and regolith, $d_{gap}\to 0$.
See section 2.1 for a detailed description of parameters used.

###  2.1 Gas pressure

The sublimation rate at the bottom of the ice slab controls both the pressure
build-up and the flow rate into the spider vent. Upon initial fracturing there
may be a very energetic flow, driven by the large initial pressure gradient,
followed by a reduction to a steady state balance between gas production and
ejection. Both these cases were simulated with a commercial hydrodynamics code
in Thomas et al. [2011a] and Thomas et al. [2011b].

The gas production rate is estimated by comparing the amount of energy
deposited by insolation at the ice-regolith boundary with the latent heat of
sublimation of CO2. Kieffer [2007] assumed a perfectly transparent ice slab
and calculate an area gas production rate of $\dot{m}=1.157\times 10^{-4}$ kg
m-2 s-1. Thomas et al. [2011b] meanwhile assume $25\%$ of the solar energy
penetrates a $z=20$ cm thick slab (based on Pilorget et al. [2011]) and arrive
at $\dot{m}=2.21\times 10^{-5}$ kg m-2 s-1.

This constant gas production during spring–time daylight hours cannot escape
and swiftly fills up all available sub-ice volume. Gas pressure will at this
point be equal to the cryostatic pressure of $P=P_{surf}+\rho_{ice}~{}gz$,
where $g=3.72$ m s-2 is the Martian gravity and $\rho_{ice}=1606$ kg m-3 is
the solid CO2 ice density. Taking $z$ as 0.7 m, as in Pilorget et al. [2011],
and $P_{surf}\approx 600$ Pa, results in sub-ice pressures of $P_{cryo}\approx
4700$ Pa. This pressure will continue building until the stress in the ice
slab above exceeds its rupture strength and it fractures. During the initial
rupture, Thomas et al. [2011a] simulated a total flow rate of $Q_{i}=0.05-0.5$
kg s-1 draining a 400 square metre area into a central vent. Subsequently,
ejection of gas reduces the pressure gradient and, in a steady state, the 400
square metre area will produce a total flow some $\sim 5$ times smaller:
$Q_{s}=0.00884$ kg s-1 for $\dot{m}=2.21\times 10^{-5}$ kg m-2 s-1.
Hydrodynamics simulations of these magnitudes of flow rate produce gas and
dust ejection velocities of $\sim 10$ ms-1, consistent with the erosion of
regolith and the spread of fine material in fans on the surface [Thomas et
al., 2011a, b].

A gas production area of 400 square metres corresponds to a circular gas
production radius of only $r_{P}=11.3$ m, somewhat smaller than the 55 m
separation (spider radius $R\sim 25$ m) between densely packed spiders found
in Hao et al. [2019]. Gas flows here therefore represent lower limits; larger
spiders will have larger mass flow rates, but the ice breaking pressure will
remain the same.

The critical breaking pressure can be estimated using the fracturing model of
Portyankina et al. [2010], which is based on an engineering model of the
bending of a thin sheet, as shown in Figure 2. Portyankina et al. [2010] give
the maximum stress at the middle of a sheet of radius $R$, and thickness $z$,
when subjected to a constant pressure $P$ over a radius of $r_{P}$. Recently,
Kaufmann et al. [2020] have measured the yield strength of CO2 ice at the
relevant temperature as $\sigma=12.3$ MPa, and so we can set the stress equal
to this value and solve for the critical gas pressure needed to cause
fracturing

$P_{crit}=\frac{8\sigma
z^{2}}{3r_{P}^{2}}\left[4-(1-\nu)\frac{r_{P}^{2}}{R^{2}}+4(1+\nu)\ln\frac{R}{r_{P}}\right]^{-1}.$
(1)

Here, $\nu_{I}=0.544$ is the Poisson’s ratio of CO2 ice [Yamashita and Kato,
1997]. It is assumed that the pressure builds up and is constant over a
length–scale $r_{P}$, comparable to the spider radius $R$ (see Fig. 2). We
favour a ratio of $r_{P}/R$ equal to one, in contrast to the value of 0.25 as
used by Portyankina et al. [2010]. We think this is more realistic because gas
pressure should build up over the whole radius of the spider, right up to its
edge where the maximum sub-ice pressure exists, and beyond-which pressure
begins to decrease again with decreasing distance towards the next spider
(Fig. 1). Nevertheless, we plot $P_{crit}$ for a range of ratios and for two
different overall spider sizes in Figure 3.

Figure 2: Model for the breaking of a thin plate as in Portyankina et al.
[2010]. Pressure $P$ builds up in a circular region of radius $r_{P}$ under a
plate of thickness $z$ and radius $R$ with simply supported edges. Figure 3:
Gas pressure required to fracture a sheet of CO2 ice of thickness $z=0.7$ m
and radius $R$ for different scales of gas buildup area $r_{P}$.

It can be seen from Figure 3 that large pressures, exceeding the cryostatic
pressure of $P_{cryo}\sim 5$ kPa, can develop before fracturing occurs. This
is especially the case for the smaller spiders (small R) and for small gas
production radii (small $r_{P}/R$), both of which can be understood as
reducing the total area that the pressure is applied to, meaning more gas
pressure must be added to achieve the same total stress in the ice. For the
densely packed spider separations observed in Hao et al. [2019] and Hao et al.
[2020] and a gas production radius of $r_{P}=0.25R\sim 6.5$ m (even smaller in
size than the case simulated by Thomas et al. [2011b] above), a total pressure
of $P_{crit}\approx 33$ kPa can build up before fracturing occurs. Thus, in
this case our maximum pressure gradient driving gas flow at the initial breach
would be $(P_{crit}-P_{surf})/r_{P}\approx 33$ kPa$/6.5$ m.

###  2.2 Permeability

The permeability of a porous medium can be expressed as (see for example Balme
and Hagermann [2006])

$K=\frac{D^{2}\phi^{3}}{72\tau(1-\phi)^{2}},$ (2)

where $D$ is the constituent particle size, $\phi$ is porosity and $\tau$ is
tortuosity. Typical values used for Martian regolith are $\phi=40\%$,
$\tau=25/12$ [Balme and Hagermann, 2006, Demidov et al., 2015], and particle
sizes between $\sim$microns and a millimetre (e.g. Morgan et al., 2018). Over
this range of particle sizes $K$ changes by several orders of magnitude,
exceeding any uncertainty in $\phi$ and $\tau$. For gasses at a low pressure
$P$ such as we are dealing with (rather than liquid flows), K must be
corrected by the ‘Klinkenberg correction’ [Klinkenberg, 1941]

$k_{G}=k+\frac{6.9}{P}k^{0.64},$ (3)

where pressure is expressed in psi and permeability in mdarcies. Conversion
back to SI units is performed using $k=cK$ and $k_{G}=cK_{G}$, with
$c=0.986923\times 10^{-3}\mu$m2.

With the above permeability, Darcy’s law for the total mass flow, $Q$, through
a cross-sectional area $A$ of porous medium, driven by a pressure gradient of
$dP/dr$, is [Ahmed, 2010]

$Q=-\frac{\rho K_{G}A}{\eta}\frac{dP}{dr},$ (4)

with the gas density $\rho$ and viscosity $\eta$ given by

$\rho=\frac{Pm_{CO2}}{k_{B}T},$ (5)

and

$\eta=\rho\lambda\sqrt{\frac{2k_{B}T}{\pi m_{CO2}}},$ (6)

respectively. Here $k_{B}$ is Boltzmann’s constant, $m_{CO2}=7.3\times
10^{-26}$ kg and $d_{CO2}=330\times 10^{-12}$ m are the molecular mass and
diameter of CO2, and the mean free path is given by

$\lambda=\frac{k_{B}T}{\sqrt{2}\pi Pd_{CO2}}.$ (7)

It is assumed that gas temperature, $T$, is equal to the sublimation
temperature at the sub-ice pressure

$T=\frac{-b}{\ln(P/100)-a},$ (8)

where we use the constants $a=23.3494$ and $b=3182.48$, as in Thomas et al.
[2011a]. Under the 0.7 m deep slab these values are $T=163$ K (around 10 K
higher than at the surface) $\rho=0.155$ kg m-3, $\lambda=9.75\times 10^{-7}$
m, and $\eta=2.12\times 10^{-5}$ Pa s.

Assuming a cylindrical geometry with a thickness of permeable regolith $d$
metres, the cross-sectional area through which gas can flow, $A$, at a
distance $r$ from the central vent, is $A=2\pi rd$. The pressure gradient,
$\frac{dP}{dr}$, can be approximated as $\Delta P/r_{P}$, where $\Delta P$ is
the pressure drop from the high-pressure, sub-ice environment to the low-
pressure vent over a length-scale $r_{P}$. Close to the central vent (i.e.
$r\ll R$), almost the entire gas flow, $Q$ from above, must pass through an
area of permeable regolith. We therefore set $r=1$ m and obtain an expression
from Eqn. 4 for the required depth of regolith in terms of this flow as

$d=\frac{Q\eta}{2\pi\rho K_{G}}\frac{r_{P}}{\Delta P}.$ (9)

Note that we have assumed a constant permeability for the regolith layer
while, in reality, porosity, and therefore permeability, is generally reckoned
to decrease with depth (see for example Morgan et al. [2018]). Decreasing
permeability will restrict gas flow rates meaning that, again, our estimate
here is a ‘strongest-case’ for flow through porous regolith.

The length–scale $r_{P}$ is the gas buildup radius of above and should be
similar in size to the spider radius $R$. Indeed, Hao et al. [2019] suggest
that the distance between adjacent spiders is determined by the pressure fall-
off, so that new fractures open (forming a new spider) at a distance from the
existing vent where pressure increases to it’s maximum sub-ice value,
$r_{P}\approx R$. In the steady state case, therefore, a further
simplification can be made. Total flow into the vent is the sum of the per-
metre gas production rate over the area of the spider $Q=\dot{m}\pi R^{2}$,
and, having made the assumption $r_{P}=R$, Eqn. 1 can be simplified and
combined with Eqn. 9 to make the dependence on spider size $R$ explicit. We
also make the additional assumption that $\Delta P\approx P_{crit}$
(neglecting the comparatively small surface pressure) so that we are using the
maximum pressure gradient possible: the ice breaking pressure, rather than the
cryostatic pressure. Then $R$ can be expressed as

$R=\left[\frac{16d\rho K_{G}\sigma z^{2}}{3\eta\dot{m}(3+\nu)}\right]^{1/5}.$
(10)

This is the maximum radius of a cylindrical volume of porous regolith which
can be effectively drained into a central vent at a steady state production
$\dot{m}$ before gas pressure at the edge exceeds ice strength and a new vent
opens.

##  3 Results and discussion

We compute the required regolith depth $d$ from Equation 9 for a range of
regolith particle sizes and two scenarios: steady state flow and just after
the initial rupture. In the steady state, we use the minimum flow rate from
above, $Q_{s}=0.00884$ kg s-1, and cryostatic pressure over a length-scale
equal to the smaller spider separations, i.e. $\Delta
P/r_{P}=(P_{cryo}-P_{surf})/r_{P}\approx 5$kPa$/25$ m. For the rupture case we
use the strongest case of $(P_{crit}-P_{surf})/r_{P}\approx 33$ kPa$/6.5$ m,
from above, as well as the slightly higher Thomas et al. [2011a] flow rate
($Q_{i}=0.05$ kg s-1).

Figure 4: Required depth of the permeable layer, $d$, to sustain the minimum
mass flow, $Q$, for various particle diameters and pressure gradients.

Figure 4 shows the results: in the steady state, a permeable regolith layer of
thickness at least $d\approx 13.4$ m is needed to support the flow rate for
the likely particle size of up to $D\approx 200~{}\mu$m [Kieffer, 2007]. This
is rather large, given the presence of ground ice and a decreasing regolith
porosity with depth. Even in the case of the strongest possible pressure
gradient just after rupture, the required thickness is still $d\approx 5.4$ m.
Only in the case of very large regolith particles ($\sim 1$ mm) and the
strongest pressure gradients does the required layer thickness fall below one
metre. An impermeable layer of water-ice bonded regolith is generally assumed
to be present in the Martian polar regions at a depth of between several
centimetres [Portyankina et al., 2010, Pilorget et al., 2011] and the seasonal
skin depth of $\sim 0.8$ m [Kieffer, 2007]. This is supported by thermal
inertia [Bandfield and Feldman, 2008] and gamma-ray and neutron spectrometry
measurements [Boynton et al., 2002, Diez et al., 2008], as well as the recent
detection of shallow-buried water ice cliffs at mid-latitudes [Dundas et al.,
2018].

Given the above, it seems unlikely that the total gas flow into a spider vent
can occur purely through a porous regolith. What then is the largest size of a
feature than can be supported this way? Figure 5 shows the results of Equation
10, where we assumed a steady state flow with the minimum gas production and a
fixed geometry, and solve for $R$ for a set of permeable regolith depths. This
is the length-scale over which a spider can effectively drain gas to a central
vent before pressure exceeds the overlying ice strength at the edge. Even in
the case of large regolith depths, the length-scales ($\sim 15$ m) are below
the observed minimum spider separation distances ($\sim 25$ m) for a particle
size below $200~{}\mu$m. For more realistic porous layer thicknesses, only
areas of a few metres in radius can be drained. Gas flow through porous
regolith therefore seems incompatible with the expected polar stratigraphy.
Even in the most amenable case of the smallest distances between araneiforms,
with the highest pressure gradients and the lowest gas flow rates, there is
insufficient depth of permeable regolith to support the flow associated with
the whole spider area venting through a single central vent. Small areas of a
few metres in radius around local vents (such as the so-called ‘baby spiders’
identified by Schwamb et al. 2018) may be drained by viscous flow through the
pores, but it seems likely that the larger (tens to hundreds of metres across)
spiders must involve the original Kieffer [2007] model of ice slab levitation
and free gas flow.

Figure 5: Maximum radius $R$, over which steady-state gas production can be
drained before pressure build–up exceeds the overlying ice strength at the
edge, for different regolith depths.

Slab levitation is likely when considering the large gas pressures, in excess
of cryostatic pressure, calculated above, as well as the small gap size needed
for relatively large flows. This required gap size $d_{gap}$ can be computed.
Assuming laminar Poiseuille flow between two parallel plates (an approximation
that likely holds away from the energetic flow at the vent itself) and a
pressure gradient of $\Delta P/R$, the mass flow rate per unit circumference
is given by $q=\frac{\rho d_{gap}^{3}}{12\eta}\frac{\Delta P}{R}$. The total
mass flow into an annulus at radius $r$ can also be expressed in terms of the
steady state area production rate as $Q(r)=\dot{m}\pi(R^{2}-r^{2})$, so that
per unit of circumference it is $\frac{\dot{m}}{2}(R^{2}/r-r)$. The two can
then be equated and solved for $d_{gap}$, the required distance between the
two plates to accommodate the flow, as

$d_{gap}=\left[\frac{6\dot{m}\eta}{\rho\Delta
P}\left(\frac{R^{3}}{r}-rR\right)\right]^{1/3}.$ (11)

Figure 6 shows the results of equation 11 for the two spider sizes from above.
The required depth is greatest for the large mass flows of the bigger spiders,
and increases towards the central vent as the gas is concentrated together in
the converging flow (it is undefined at $r=0$, the vent itself). Nonetheless,
the required gap between the regolith and ice, $d_{gap}$ is much smaller than
the depth of permeable regolith $d$ calculated above and rarely exceeds one
centimetre. Thus the whole flow can be supported by a gap of a few
centimetres, as suggested by Kieffer [2007].

Figure 6: Required gap size, $d_{gap}$, to support free-flow of the total gas
production, between the ice and regolith, across the radii of two different
spider sizes.

Such a regolith–ice gap may occur due to levitation of the ice slab, but also
important is the fact that the slab will bend upwards before fracturing,
providing an additional gap for gas to flow through. Portyankina et al. [2010]
give an expression for the maximum displacement at the centre of the ice sheet
(see Fig. 2), which we can evaluate just before rupture by using the pressure
$P_{crit}$ calculated above:

$w_{crit}=\frac{3(1-\nu)r_{P}^{2}P_{crit}}{16Ez^{3}}\left[4(3+\nu)R^{2}-(7+3\nu)r_{P}^{2}-4(1+\nu)r_{P}^{2}\ln\frac{R}{r_{P}}\right],$
(12)

with the CO2 ice Young’s modulus of $E=11.5$ GPa [Yamashita and Kato, 1997].
The resulting displacement is shown in Figure 7, for the same spider sizes and
$r_{P}/R$ ratios as before. The ice in the centre of 25 m and 150 m radius
spiders is bent upwards by $\sim 20$ cm and $\sim 10$ m, respectively, before
fracturing occurs. These values greatly exceed the required gap sizes for the
gas flow shown in Fig. 6. Therefore, even though the upwards displacement is
smaller towards the edge of the feature, it seems likely that upwards bending
of the ice sheet can contribute significantly to the opening of a regolith–ice
gap, aiding the flow of gas towards the spider vent.

Figure 7: Maximum upwards displacement in the centre of an ice sheet just
before fracturing, for different sizes $R$, and gas buildup areas $r_{P}$.

With viscous flow through regolith pores unable to deliver the full gas
production of a reasonable sized spider, and large gas pressure enough to bend
or levitate the ice slab, it may be that some spiders undergo both types of
flow. For example, low volumes of gas moving through a large cross-sectional
area of porous regolith at the outskirts of a feature or around multiple
vents, transitioning to free-flow nearer the vent where the converging gas has
bent or levitated the overlying ice slab upwards, opening up a gap. Meanwhile,
flow through and just above a loose regolith will certainly begin to erode
material, digging troughs downwards into it. In Hao et al. [2019]’s conceptual
model, regolith properties are cited as explaining the different appearances
and spacing of spider sub-types. Although the above calculations demonstrate
that significant gas flow through the pores is unlikely, regolith properties,
such as cohesion, could still have an important effect on morphology. So-
called ‘fat spiders’ could still be examples of low-cohesion, more easily
eroded regoliths, as suggested by Hao et al. [2019], for example. Once these
troughs and pits are eroded into the regolith, gas will preferentially flow
down and further erode them, as described in Hao et al. [2019]. This will
provide additional flow paths for the gas but, since slab ice is thought to
conformally cover the terrain in winter, at least some levitation or bending
is required to open up a gap again in the spring. The physics of mixed gas
flow through the regolith pores, eroded troughs, and a levitated or bent-open
ice gap will be complicated, and may involve turbulent as well as viscous
flow. Alternatively, local variations in the depth of the impermeable water-
ice layer, below the resolution of orbital observations, may allow increased
porous gas flow in some localised areas. It seems unlikely, however, that
water-ice depths could vary from $\sim$cm to $\sim$metres over a few metres
lateral distance in otherwsie identical terrain.

Finally, it should also be noted that some recent work [Chinnery et al., 2018]
suggests a much lower translucence of CO2 ice when compared to previous
experiments [Hansen, 1999], which would require much thinner ice sheets (on
the order of 10 cm) in order to trigger the solid state greenhouse effect. We
note that in this case, of a reduced $z$, cryostatic pressure and the critical
breaking pressure would be reduced (see Eqn. 1), further lowering the pressure
gradient and making permeable gas flow even less effective, even for small gas
flow rates.

##  4 Conclusion

In the formation of the so–called araneiform features, seen in Mars’ southern
polar regions, subsurface flow of CO2 gas below an ice slab can occur in a gap
between ice and regolith, and/or in the substrate itself. The latter
possibility has recently received increased attention because spider spatial
density scales seem to be related to properties of the underlying regolith. We
investigate the effectiveness of subsurface flow of CO2 and compare it to gas
flow through the gap between ice and regolith. Based on previously estimated
flow rates and thermophysical arguments, we suggest that there is insufficient
depth of porous regolith to support the full gas flow of all but the very
smallest observed spiders through the regolith. On the other hand, free gas
flow through a regolith–ice gap is capable of supplying the likely flow rates
for gap sizes on the order of a centimetre. This size of gap can be opened in
the centre of a spider feature by gas pressure bending the overlying ice slab
upwards, or by levitating it entirely as suggested in the original Kieffer
[2007] model. Our calculations therefore support at least some of the gas
flowing through a gap opened between the regolith and ice. Regolith properties
most likely still play a role in the evolution of spider morphology, by
regolith cohesion controlling the erosion of the central pit and troughs, for
example.

## Acknowledgements

The authors would like to thank Jingyan Hao and an anonymous reviewer for
their constructive and helpful comments which helped improve this manuscript
substantially. NA and AH acknowledge support from the UK Space Agency, grant
no. ST/R001375/2. EK was supported by STFC grant no. ST/S001271/1.

## References

  * Ahmed [2010] Ahmed, T., 2010.  Reservoir Engineering Handbook.  Fourth ed., Gulf Professional Publishing, Boston.
  * Balme and Hagermann [2006] Balme, M., Hagermann, A., 2006\.  Particle lifting at the soil-air interface by atmospheric pressure excursions in dust devils.  Geophysical Research Letters 33, L19S01.  doi:10.1029/2006GL026819.
  * Bandfield and Feldman [2008] Bandfield, J., Feldman, W., 2008\.  Martian high latitude permafrost depth and surface cover thermal inertia distributions.  Journal of Geophysical Research: Planets 113, 0148–0227.
  * Boynton et al. [2002] Boynton, W.V., Feldman, W.C., Squyres, S.W., Prettyman, T.H., Brückner, J., Evans, L.G., Reedy, R.C., Starr, R., Arnold, J.R., Drake, D.M., Englert, P.A.J., Metzger, A.E., Mitrofanov, I., Trombka, J.I., d’Uston, C., Wänke, H., Gasnault, O., Hamara, D.K., Janes, D.M., Marcialis, R.L., Maurice, S., Mikheeva, I., Taylor, G.J., Tokar, R., Shinohara, C., 2002.  Distribution of hydrogen in the near surface of mars: Evidence for subsurface ice deposits.  Science 297, 81–85.  doi:10.1126/science.1073722.
  * Brown and Matson [1987] Brown, R.H., Matson, D.L., 1987\.  Thermal effects of insolation propagation into the regoliths of airless bodies.  Icarus 72, 84–94.  doi:10.1016/0019-1035(87)90122-9.
  * Chinnery et al. [2018] Chinnery, H.E., Hagermann, A., Kaufmann, E., Lewis, S.R., 2018\.  The Penetration of Solar Radiation Into Carbon Dioxide Ice.  Journal of Geophysical Research (Planets) 123, 864–871.  doi:10.1002/2018JE005539.
  * Demidov et al. [2015] Demidov, N.E., Bazilevskii, A.T., Kuz’min, R.O., 2015.  Martian soils: Varieties, structure, composition, physical properties, drillability, and risks for landers.  Solar System Research 49, 209–225.  doi:10.1134/S0038094615040024.
  * Diez et al. [2008] Diez, B., Feldman, W., Maurice, S., Gasnault, O., Prettyman, T., Mellon, M., Aharonson, O., Schorghofer, N., 2008\.  H layering in the top meter of mars.  Icarus 196, 409 – 421.  doi:https://doi.org/10.1016/j.icarus.2008.02.006.
  * Dundas et al. [2018] Dundas, C.M., Bramson, A.M., Ojha, L., Wray, J.J., Mellon, M.T., Byrne, S., McEwen, A.S., Putzig, N.E., Viola, D., Sutton, S., Clark, E., Holt, J.W., 2018\.  Exposed subsurface ice sheets in the martian mid-latitudes.  Science 359, 199–201.  doi:10.1126/science.aao1619.
  * Fanale et al. [1990] Fanale, F., Salvail, J., Matson, D., Brown, R., 1990\.  The effect of volume phase changes, mass transport, sunlight penetration, and densification on the thermal regime of icy regoliths.  Icarus 88, 193–204.  doi:10.1016/0019-1035(90)90185-C.
  * Hansen et al. [2010] Hansen, C.J., Thomas, N., Portyankina, G., McEwen, A., Becker, T., Byrne, S., Herkenhoff, K., Kieffer, H., Mellon, M., 2010.  HiRISE observations of gas sublimation-driven activity in Mars’ southern polar regions: I. Erosion of the surface.  Icarus 205, 283–295.  doi:10.1016/j.icarus.2009.07.021.
  * Hansen [1999] Hansen, G.B., 1999.  Control of the radiative behavior of the Martian polar caps by surface CO2 ice: Evidence from Mars Global Surveyor measurements.  JGR Planets 104, 16471–16486.  doi:10.1029/1998JE000626.
  * Hansen [2005] Hansen, G.B., 2005.  Ultraviolet to near-infrared absorption spectrum of carbon dioxide ice from 0.174 to 1.8 $\mu$m.  Journal of Geophysical Research (Planets) 110, E11003.  doi:10.1029/2005JE002531.
  * Hao et al. [2019] Hao, J., Michael, G.G., Adeli, S., Jaumann, R., 2019\.  Araneiform terrain formation in Angustus Labyrinthus, Mars.  Icarus 317, 479–490.  doi:10.1016/j.icarus.2018.07.026.
  * Hao et al. [2020] Hao, J., Michael, G.G., Adeli, S., Jaumann, R., Portyankina, G., Hauber, E., Millot, C., Zuschneid, W., 2020\.  Variability of spider spatial configuration at the Martian south pole.  Planet. Space Sci. 185, 104848\.  doi:10.1016/j.pss.2020.104848.
  * Kaufmann et al. [2020] Kaufmann, E., Attree, N., Bradwell, T., Hagermann, A., 2020\.  Hardness and Yield Strength of CO2 Ice Under Martian Temperature Conditions.  JGR Planets 125, e06217.  doi:10.1029/2019JE006217.
  * Kieffer [2000] Kieffer, H.H., 2000.  Annual Punctuated CO2 Slab-Ice and Jets on Mars, in: Second International Conference on Mars Polar Science and Exploration, p. 93.
  * Kieffer [2007] Kieffer, H.H., 2007.  Cold jets in the Martian polar caps.  Journal of Geophysical Research (Planets) 112, E08005.  doi:10.1029/2006JE002816.
  * Kieffer et al. [2006] Kieffer, H.H., Christensen, P.R., Titus, T.N., 2006.  CO2 jets formed by sublimation beneath translucent slab ice in Mars’ seasonal south polar ice cap.  Nature 442, 793–796.  doi:10.1038/nature04945.
  * Klinkenberg [1941] Klinkenberg, L.J., 1941.  The permeability of porous media to liquids and gases.  Drilling and Production Practice .
  * Martínez et al. [2012] Martínez, G.M., Renno, N.O., Elliott, H.M., 2012.  The evolution of the albedo of dark spots observed on Mars polar region.  Icarus 221, 816–830.  doi:10.1016/j.icarus.2012.09.008.
  * Mellon et al. [2004] Mellon, M.T., Feldman, W.C., Prettyman, T.H., 2004.  The presence and stability of ground ice in the southern hemisphere of Mars.  Icarus 169, 324–340.  doi:10.1016/j.icarus.2003.10.022.
  * Morgan et al. [2018] Morgan, P., Grott, M., Knapmeyer-Endrun, B., Golombek, M., Delage, P., Lognonné, P., Piqueux, S., Daubar, I., Murdoch, N., Charalambous, C., Pike, W.T., Müller, N., Hagermann, A., Siegler, M., Lichtenheldt, R., Teanby, N., Kedar, S., 2018.  A Pre-Landing Assessment of Regolith Properties at the InSight Landing Site.  Space Science Reviews 214, 104\.  doi:10.1007/s11214-018-0537-y.
  * Pilorget and Forget [2016] Pilorget, C., Forget, F., 2016\.  Formation of gullies on Mars by debris flows triggered by CO2 sublimation.  Nature Geoscience 9, 65–69.  doi:10.1038/ngeo2619.
  * Pilorget et al. [2011] Pilorget, C., Forget, F., Millour, E., Vincendon, M., Madeleine, J.B., 2011.  Dark spots and cold jets in the polar regions of Mars: New clues from a thermal model of surface CO 2 ice.  Icarus 213, 131–149.  doi:10.1016/j.icarus.2011.01.031.
  * Piqueux et al. [2003] Piqueux, S., Bryne, S., Richardson, M.I., 2003.  Sublimation of Mars’s southern seasonal CO2 ice cap and the formation of spiders.  JGR Planets 104, 5084\.  doi:10.1029/2002JE002007.
  * Piqueux and Christensen [2008] Piqueux, S., Christensen, P.R., 2008\.  North and south subice gas flow and venting of the seasonal caps of Mars: A major geomorphological agent.  JGR Planets 113, E06005.  doi:10.1029/2007JE003009.
  * Portyankina et al. [2010] Portyankina, G., Markiewicz, W.J., Thomas, N., Hansen, C.J., Milazzo, M., 2010.  HiRISE observations of gas sublimation-driven activity in Mars’ southern polar regions: III. Models of processes involving translucent ice.  Icarus 205, 311–320.  doi:10.1016/j.icarus.2009.08.029.
  * Schwamb et al. [2018] Schwamb, M.E., Aye, K.M., Portyankina, G., Hansen, C.J., Allen, C., Allen, S., Calef, F.J., Duca, S., McMaster, A., Miller, G.R.M., 2018\.  Planet Four: Terrains - Discovery of araneiforms outside of the South Polar layered deposits.  Icarus 308, 148–187.  doi:10.1016/j.icarus.2017.06.017, arXiv:1708.07858.
  * Tanaka et al. [2014] Tanaka, K.L., Skinner, J. A., J., Dohm, J.M., Irwin, R. P., I., Kolb, E.J., Fortezzo, C.M., Platz, T., Michael, G.G., Hare, T.M., 2014.  Geologic map of Mars.  U.S. Geological Survey Report 3292\.  doi:10.3133/sim3292.
  * Thomas et al. [2010] Thomas, N., Hansen, C.J., Portyankina, G., Russell, P.S., 2010\.  HiRISE observations of gas sublimation-driven activity in Mars’ southern polar regions: II. Surficial deposits and their origins.  Icarus 205, 296--310.  doi:10.1016/j.icarus.2009.05.030.
  * Thomas et al. [2011a] Thomas, N., Portyankina, G., Hansen, C.J., Pommerol, A., 2011a.  HiRISE observations of gas sublimation-driven activity in Mars’ southern polar regions: IV. Fluid dynamics models of CO2 jets.  Icarus 212, 66--85.  doi:10.1016/j.icarus.2010.12.016.
  * Thomas et al. [2011b] Thomas, N., Portyankina, G., Hansen, C.J., Pommerol, A., 2011b.  Sub-surface CO2 gas flow in Mars’ polar regions: Gas transport under constant production rate conditions.  Geophysical Research Letters 38, L08203.  doi:10.1029/2011GL046797.
  * Yamashita and Kato [1997] Yamashita, Y., Kato, M., 1997\.  Viscoelastic properties of polycrystalline solid methane and carbon dioxide.  Geophys. Res. Lett. 24, 1327--1330.  doi:10.1029/97GL01205.