General scaling relations for locomotion in granular media
James Slonaker$^1$
D. Carrington Motley$^1$
Qiong Zhang$^1$
Stephen Townsend$^1$
Carmine Senatore$^1$
Karl Iagnemma$^1$
Ken Kamrin
%{Using Resistive Force Theory, a reduced-order empirical model for granular drag, we have discovered fundamental scaling laws, which instruct how to scale size, mass, and driving parameters to model forces generated during locomotion. We confirm these relations by testing pairs of 3D-printed wheels of various shapes and sizes under many load and motion conditions in a controlled sand bed. We also provide an explanation for how the scalings may be obtained through the assumptions of Coulomb plasticity, and how the relations can be augmented to predict wheel performance under different gravities. %{Recent experiments have demonstrated the validity of Resistive Force Theory (RFT), a reduced-order empirical model for granular drag, to predict dynamical interactions of solid bodies with granular materials. Through a derivation utilizing certain RFT hypotheses, we have discovered the first family of fundamental scaling laws in granular locomotion, which instruct how to down- or up-scale size, mass, and driving parameters in order to scale-model forces generated during the locomotion process. The relative simplicity of RFT asserts that only one material parameter, the ``grain-structure coefficient'', is required, which reduces the complexity of the non-dimensionalization of the system considerably. We experimentally confirm these relations by testing pairs of 3D-printed wheels of various shapes and sizes under many load and motion conditions in a controlled sand bed. We also provide an explanation for how the scalings may be obtained alternatively through the assumptions of Coulomb plasticity, and how the relations can be further augmented to predict wheel performance under different ambient gravity. {We derive a general dimensionless form for granular locomotion, which is validated in experiments and Discrete Element Method (DEM) simulations. The form instructs how to scale size, mass, and driving parameters in order to relate dynamic behaviors of different locomotors in the same granular media. The scaling can be derived by assuming intrusion forces arise from Resistive Force Theory (RFT) or equivalently by assuming the granular material behaves as a continuum obeying a frictional yield criterion. The scalings are experimentally confirmed using pairs of wheels of various shapes and sizes under many driving conditions in a common sand bed. We discuss why the two models provide such a robust set of scaling laws even though they neglect a number of the complexities of granular rheology. Motivated by potential extra-planetary applications, the dimensionless form also implies a way to predict wheel performance in one ambient gravity based on tests in a different ambient gravity. We confirm this using DEM simulations, which show that scaling relations are satisfied over an array of driving modes even when gravity differs between scaled tests.
RFT_ReSubmit3.tex
\documentclass[twocolumn,english,aps,prl,superscriptaddress]{revtex4}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{snapshot}
\makeatletter

\usepackage{amsfonts}
\usepackage{amsthm}
\usepackage{mathrsfs}
\usepackage{graphics}
\usepackage{lscape}
\usepackage{changebar}
\usepackage{subfigure}
%\usepackage{subcaption}
\usepackage{epsfig}
\usepackage{babel}
\usepackage{array}
\usepackage{color}
\usepackage{cmap}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{mathtools}
\usepackage[T1]{fontenc}
\newcommand{\edit}[1]{{\color{black}{#1}}}
%\newcommand{\jcom}[1]{{\color{red}{[KEN: #1]}}}
%\newcommand{\check}[1]{{\color{red}{[#1]}}}

\bibliographystyle{apsrev}

\makeatother

\begin{document}

\title{General scaling relations for locomotion in granular media}

\author{James Slonaker$^1$, D. Carrington Motley$^1$, Qiong Zhang$^1$, Stephen Townsend$^1$, Carmine Senatore$^1$, Karl Iagnemma$^1$, Ken Kamrin}
 \email{kkamrin@mit.edu}
 \affiliation{Department of Mechanical Engineering, MIT, Cambridge, MA 02139, USA}


\date{\today}

\begin{abstract}
%{Using Resistive Force Theory, a reduced-order empirical model for granular drag, we have discovered fundamental scaling laws, which instruct how to scale size, mass, and driving parameters to model forces generated during locomotion. We confirm these relations by testing pairs of 3D-printed wheels of various shapes and sizes under many load and motion conditions in a controlled sand bed. We also provide an explanation for how the scalings may be obtained through the assumptions of Coulomb plasticity, and how the relations can be augmented to predict wheel performance under different gravities.}
%{Recent experiments have demonstrated the validity of Resistive Force Theory (RFT), a reduced-order empirical model for granular drag, to predict dynamical interactions of solid bodies with granular materials.  Through a derivation utilizing certain RFT hypotheses, we have discovered the first family of fundamental scaling laws in granular locomotion, which instruct how to down- or up-scale size, mass, and driving parameters in order to scale-model forces generated during the locomotion process.   The relative simplicity of RFT asserts that only one material parameter, the ``grain-structure coefficient'', is required, which reduces the complexity of the non-dimensionalization of the system considerably.  We experimentally confirm these relations by testing pairs of 3D-printed wheels of various shapes and sizes under many load and motion conditions in a controlled sand bed. We also provide an explanation for how the scalings may be obtained alternatively through the assumptions of Coulomb plasticity, and how the relations can be further augmented to predict wheel performance under different ambient gravity.}

{We derive a general dimensionless form for granular locomotion, which is validated in experiments and Discrete Element Method (DEM) simulations.  The form instructs how to scale size, mass, and driving parameters in order to relate dynamic behaviors of different locomotors in the same granular media.  The scaling can be derived by assuming intrusion forces arise from Resistive Force Theory (RFT) or equivalently by assuming the granular material behaves as a continuum obeying a frictional yield criterion.   The scalings are experimentally confirmed using pairs of wheels of various shapes and sizes under many driving conditions in a common sand bed. We discuss why the two models provide such a robust set of scaling laws even though they neglect a number of the complexities of granular rheology.  Motivated by potential extra-planetary applications, the dimensionless form also implies a way to predict wheel performance in one ambient gravity based on tests in a different ambient gravity. We confirm this using DEM simulations, which show that scaling relations are satisfied over an array of driving modes even when gravity differs between scaled tests.}


\end{abstract}

\pacs{PACS number(s): }

\maketitle

\section{Introduction}
Due to the complexity of the constitutive behavior of granular media \cite{goddard2014continuum} dynamic interactions between grains and solid bodies are challenging to model without resorting to grain-by-grain discrete particle methods.  When materials have a simple constitutive behavior, one benefit is the ability to develop straightforward scaling laws, which can be used to study material dynamics in complicated geometries. Newtonian fluids, for example, have a well-known rheology governed by two material parameters, viscosity and density, which leads to a handful of well-known dimensionless parameters (e.g. Reynolds number, Froude number) that can be used for controlled scaling analyses.  In contrast, granular media are more complex and the constitutive behavior is still widely debated.  

Granular materials display a diversity of phenomena that require special constitutive treatment to model, such as history- and preparation-dependent strengthening and dilation \cite{schofield1968critical,wood1990soil,roux1998texture}, flow anisotropy and normal stress differences \cite{depken07,weinhart2013coarse,luding08}, nonlinear rate-sensitive yielding \cite{bagnold54,dacruz05,midi2004dense,jop06}, and nonlocality due to the finite size of grains \cite{kamrin2012nonlocal,henann2013predictive,kamrin2015nonlocal,mohan02,aranson02}. Obtaining pragmatic scaling relations from these models or a combination thereof has some inherent difficulties: (1) A quantitative model that attempts to capture the various granular complexities invariably requires more material parameters. These invoke more dimensionless numbers, each of which acts as a constraint on how to make a scaled pair of experiments. (2) Many granular continuum models reference grain properties directly, such as the mean grain diameter. This can lead to scaling relations that require scaled tests to be performed with different grains.  The ability to manufacture such granular systems in practice would be highly nontrivial, and moreover, changing grain size can introduce new physics to the grain-grain interaction, such as the relative importance of charge as grains enter the powder regime. 

Despite the complexities of the constitutive behavior of grains, recent results suggest that the specific problem of calculating resistive forces on generally-shaped, rigid intruders is often well-described with a simple model called the  Resistive Force Theory (RFT) \cite{li13}.  RFT is an empirical model utilizing a set of hypotheses about local drag forces to approximate resistance on general solid surfaces moving in granular soils near the surface.  RFT was initially developed for viscous drag problems \cite{lighthill75}, however it has shown surprising effectiveness in granular media, where it has been used to simulate the dynamics of legged reptiles and robots \cite{li13}, swimming sandfish \cite{goldman14}, and the distribution of lift forces on curved bodies submerged in grains \cite{ding11}.  

In light of its effectiveness in multiple geometries and its dependence on very few model parameters, a natural question to ask is whether granular RFT, when combined with Newton's laws, produces a set of intruder dynamics possessing \edit{scaling behaviors}.  If these could be identified and validated experimentally, they would provide a physical basis to directly relate different granular locomotion problems in the same soil without the need to run any models or perform any discrete simulation. In application, they could be exploited as scaling laws to predict the performance of a ``large'' locomotor in a sand bed --- such as a truck wheel or a tank --- by appropriately down-scaled analysis of a smaller locomotor in the same bed.  Such a capability could be leveraged to produce scaling methods for soil interactions similar to those exploited in aerodynamics and hydrodynamics.  

Herein, we study arbitrarily shaped wheels and derive and experimentally validate a family of geometrically-general scaling laws governing their driving behaviors.  The relations are obtained through analysis of the RFT terradynamical system. As a secondary justification, we show that the same scalings arise by modeling the grains as a simple, rate-independent, frictional continuum.  This raises important questions about the importance of the constitutive complexities of granular flow when studying the specific problem of granular intrusion forces.  It also supports recent results \cite{askari2016intrusion} showing that RFT modeling assumptions emerge strongly in intrusion through a hypothetical frictional-plastic media. A key point about the scaling relations we propose is that they do not require changing the grains themselves in order to relate two experiments with different locomotors. 

The analysis we provide could be applied to locomotors with more moving parts, but here we test the concept on solid wheels for simplicity.  The scaling forms admit variation in the ambient gravity as well, which could have a number of applications in extra-planetary studies; this would admit one to perform a test on earth that predicts how a locomotor functions, say, on Mars.  This ability is validated with discrete element method (DEM) tests, which enable control over gravity.  They show that the scaling relations work over a wide range of locomotive behaviors even as gravity varies between pairs of scaled tests. 

\begin{figure}[t!!!!!!]
\includegraphics[width=0.74\columnwidth,trim={5.8cm 8.6cm 6.4cm 7.6cm},clip]{Fig1CombTry2.pdf}
%%\includegraphics[width=0.96\columnwidth,trim={3.9cm 11.6cm 3.6cm 11.6cm},clip]{Fig1Comb4.pdf}
%%\includegraphics[width=0.2\textwidth,trim={0cm 0cm 1.8cm 0cm},clip]{GeneralWheel4.pdf}\ 
%%\includegraphics[width=0.225\columnwidth,trim={0.cm 0.55cm 0.0cm 0.cm},clip]{CylDef2.pdf}
%%\includegraphics[width=0.3\columnwidth,trim={0cm 0cm 0.0cm 0cm},clip]{LugDef.pdf}%
%\begin{center}
%\includegraphics[width=0.27\textwidth,trim={0cm 0cm 1.8cm 0cm},clip]{GeneralWheel4.pdf}\ \end{center}
%%\newline
%\begin{center}
%\includegraphics[width=0.295\columnwidth,trim={0.cm 0.55cm 0.0cm 0.cm},clip]{CylDef2.pdf}
%\
%\includegraphics[width=0.3925\columnwidth,trim={0cm 0cm 0.0cm 0cm},clip]{LugDef.pdf}
%\end{center}%
%%\caption{\edit{(a) General problem: Driving an arbitrarily shaped wheel of width $D$ and rotational velocity $\omega$, carrying a mass $M$ under gravity $g$.  The surface element's angle of tilt $\beta$ and angle of velocity $\gamma$ are defined from the horizontal.} (b) The two wheel shapes, $f$, used in our study: cylinders (left) and lugs (right), and their corresponding definitions of $L$.}
\caption{\edit{(a) General problem: Driving an arbitrarily shaped wheel of width $D$ and rotational velocity $\omega$, carrying a mass $M$ under gravity $g$. Blow-up of a surface element shown. (b),(c) The two wheel shapes used in our study and their corresponding definitions of $L$.}}
\label{genwheel}
\end{figure}

\section{Resistive Force Theory}
Granular RFT is concerned solely with computing the resistive force on a body moving through a gravitationally loaded bed of grains.  It does not attempt to describe the flow or stress fields in the granular media.  Its basic premise, which is not actually derived physically, though new efforts are being made to provide its foundations \cite{askari2016intrusion}, is that \edit{the stress $(\sigma_x,\sigma_z)$ acting on any surface element of an intruder is determined only by the element's motion direction ($\gamma$), tilt ($\beta$), and depth $|z|$ (see Fig \ref{genwheel}(a)).  Particularly, it is assumed that $\sigma_x=\alpha_x(\beta,\gamma)\cdot\xi \cdot|z|$ and $\sigma_z=\alpha_z(\beta,\gamma)\cdot\xi \cdot|z|$ where $\alpha_x$ and $\alpha_z$ are dimensionless empirical functions with a universal known shape across many materials \cite{li13}, and $\xi$ is a ``grain-structure coefficient.''}  With units of force per volume, $\xi$ is the only parameter needed to execute an RFT calculation.  It depends on the grains, the intruder roughness, and the value of gravity.  By contrast, common engineering terradynamics models like that of Wong and Reece (based on Bekker's work) \cite{wong67} or the NATO Reference Mobility Model \cite{ahlvin92} require a larger number of fit parameters, though they model more than resistive forces.  Unlike Bekker's model \cite{bekker}, no knowledge of the material deformation is required for RFT.
%In a 3D sand bed, assume a quasi-2D intruder (e.g. a driven lug wheel, per Fig \ref{genwheel}(a)) whose motion is in the $xz$-plane, where gravity points in $-\hat{z}$ direction and $z=0$ represents the free surface. The resistive force $\mathbf{F}^{\text{res}}=(F^{\text{res}}_x,F^{\text{res}}_z)$ acting on the leading surface of the intruder, $S$, is modeled under RFT to obey
%\edit{\begin{equation}\label{RFTeqn}
%(F^{\text{res}}_x,F^{\text{res}}_z)=\int_S \Big(\alpha_x(\beta,\gamma),\ \alpha_z(\beta,\gamma)\Big) \ H(z)\ |z| \ dA_s
%\end{equation}}
%where $dA_s$ is the area of a surface element, $\beta$ is the surface element's angle of tilt, and $\gamma$ is the angle of the velocity of the surface element. Both $\beta$ and $\gamma$ are measured from the horizontal. The Heaviside function, $H(z)$, removes resistive force on parts of the intruder outside the bed, and the dependence on $|z|$ models added resistance with depth due to gravity.  \edit{The above essentially states that the force on each surface element is given by its motion direction ($\gamma$), tilt ($\beta$), and depth (see Fig \ref{genwheel}(a))}. 
%The key constitutive ingredient in the model is the selection of the functions $\alpha_x$ and $\alpha_z$, which is done empirically using experimental data on small intruding flat plates under various $\beta$ and $\gamma$ conditions. 
%Also the assumption of linear variation of surface stress with $z$; during intrusion, significant variation from a lithostatic pressure distribution is observed \cite{brzinski2013depth,guillard2014lift}.  
%
%\begin{figure}
%\includegraphics[width=0.35\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{AngleDef.pdf}
%%(a) \includegraphics[width=0.2\textwidth,trim={0cm 1.5cm 0.4cm 0cm},clip]{GeneralWheel2.pdf}\ (b)\includegraphics[width=0.45\columnwidth,trim={0cm 0cm 0cm 0cm},clip]{Cyl4Bar.pdf}
%\caption{\edit{The surface element's angle of tilt $\beta$ and angle of velocity $\gamma$ are measured from the horizontal.}}
%\label{angledef}
%\end{figure}
%It was observed \cite{li13} that the \edit{empirical functions at the heart of RFT} %$\alpha_{x}$ and $\alpha_{z}$ functions 
%have a fairly universal shape across many materials, and hence they can be determined from a single scaling factor, $\xi$, denoted the ``grain-structure coefficient.''  With units of force per volume, $\xi$ is the only parameter needed to execute an RFT calculation. It depends on the grains, the intruder roughness, and the value of gravity.  By contrast, common engineering terradynamics models like that of Wong and Reece (based on Bekker's work) \cite{wong67} or the NATO Reference Mobility Model \cite{ahlvin92} require a much larger number of fit parameters, though they model more than resistive forces.  Unlike Bekker's model \cite{bekker}, no knowledge of the material deformation is required for RFT.
%
%Resistive Force Theory is based on a set of experimental measurements of the intrusion force acting on a small, submerged, flat plate as it moves within a gravitationally loaded granular bed.  In The plate is given various angles of attack $\beta$ (i.e. the direction the plate moves) and angles of intrusion $\gamma$ (i.e. the angle the plate makes with the free surface) and the results are compiled into an empirical .  Through comparing the experimental data for $\alpha$ a generic fitting function $ \alpha_{z,x}^\text{gen}$ was found that is scalable to other non-cohesive granular materials by a constant, $\xi$, known as the ``grain-structure coefficient,'' whose value depends on the surface roughness of the solid and the properties of the granular media.  To apply RFT to simulate the forces acting on an arbitrary shaped intruder, the intruder is first modeled as a collection of discrete planar bodies.  The total resistive force $F_{z,x}$ acting on each body $s$ is then approximated as
%\begin{equation}
%F_{z,x}=\int_s \xi \alpha_{z,x}^\text{gen} \left(\beta_s, \gamma_s\right) |z|_s dA_s
%\label{RFTeqn}
%\end{equation}
%where $|z|$ is the depth below the surface and $dA$ is the area \cite{li13}.  Utilizing the linear superposition principle, the total resistive force on the intruder is approximated as the sum of the resistive forces on each planar body.  Within RFT, the only material parameter is the scaling constant or grain structure coefficient $\xi$, which is experimentally determined for each material and used to scale the fitting function.  By contrast, the common terradynamics models in use today, such as the Wong model \cite{wong67} or the Nato Reference Mobility Model \cite{ahlvin92}, all require a much larger number of fit parameters.  Due to this simplicity, we believed dimensional analysis could be used to derive fundamental scaling relations previously unknown.

\section{Dimensional Analysis}
%\noindent{\textit{Dimensional analysis:}} 
The generic wheel has a dimensionless shape, which we denote by a point-set $f$, a constant width $D$ into the plane, a characteristic length $L$ that scales the shape $f$ to give the actual wheel cross-section, and a mass $M$ assumed concentrated on the axle.  The wheel is given a fixed rotational velocity $\omega$, is acted upon by some gravity $g$, and interacts with the sand bed through some grain-structure coefficient $\xi$.  The outputs we are interested in are the power expended as the wheel drives in granular media, $P$, and the wheel's  $x$-translational velocity $V$, although other outputs could be considered.  Both outputs are functions of time, $t$. Before applying dimensional analysis, note that %Eq \ref{RFTeqn} implies the problem's dependence on $\xi$ and $D$ is seen only through the product $\xi D$.  
\edit{since the wheel surfaces do not vary in the through-thickness direction, when integrating $(\sigma_x,\sigma_z)$ over the wheel surface the dependence of the resistive force on $\xi$ and $D$ is seen only through the product $\xi D$.}  With this and a standard nondimensionalization, the wheel's steady driving limit-cycle is predicted to obey the form:
\begin{equation}
%\left[\frac{P}{mL^2\omega^3},\frac{V}{L\omega}\right]=\boldsymbol{\Psi}_\text{RFT}\left(\omega t,f,\frac{g}{L\omega^2},\frac{\xi D L^2}{mg}\right)
\left[\frac{P}{Mg\sqrt{Lg}},\frac{V}{\sqrt{Lg}}\right]=\boldsymbol{\Psi}\left(\sqrt{\frac{g}{L}} t,f,\frac{g}{L\omega^2},\frac{\xi D L^2}{Mg}\right)
\label{RFTscalND}
\end{equation}
$\boldsymbol{\Psi}$ is a four-input, two-output function as shown. If the gravity, wheel surface texture, and granular media are fixed, $g$ and $\xi$ can be absorbed into the function, giving the reduced form:
\begin{equation}
%\left[\frac{P}{mL^2\omega^3},\frac{V}{L\omega}\right]=\boldsymbol{\Psi}_\text{RFT}\left(\omega t,f,\frac{1}{L\omega^2},\frac{D L^2}{m}\right)
\left[\frac{P}{M\sqrt{L}},\frac{V}{\sqrt{L}}\right]=\tilde{\boldsymbol{\Psi}}\left(\sqrt{\frac{1}{L}} t,f,\frac{1}{L\omega^2},\frac{D L^2}{M}\right)
\label{RFTscalD}
\end{equation}
The above implies the following family of scaling relations: Consider two experiments with the same $f$, $g$, and $\xi$, but one has inputs $(L,M, D,\omega)$  and the other has inputs $( L',M',D',\omega')=( r L, s M, sr^{-2}D, r^{-1/2}\omega)$ for any positive scalars $r$ and $s$. Then the corresponding driving cycles should obey $\langle P'\rangle=s r^{1/2}\langle P\rangle$ and $\langle V'\rangle=r^{1/2}\langle V\rangle$  where $\langle\cdot\rangle$ denotes a time-average.\\

\section{Connection to Coulomb Plasticity}
%\noindent{{\textit{Connection to Coulomb plasticity:}
Although RFT is empirical, interestingly, Eq \ref{RFTscalD} can also be deduced mechanically by considering wheel motion in a 3D bed of ideal Coulomb material \cite{nedderman}. Here, the grains are treated as a rate-independent frictionally yielding continuum. Such a model could be used to predict the entire sand motion field, but instead consider dimensional analysis implications that can be identified without solving for flow. The inputs to the problem are the wheel descriptors and now three material parameters from the continuum model: the density of the material $\rho$, the material's coefficient of internal friction $\mu$, and the coefficient of sliding friction of the wheel-material interface $\mu_w$.  Hence, wheels driving through this hypothetical continuum must obey the following dimensionless form:
\begin{equation}
%\left[\frac{P}{mL^2\omega^3},\frac{V}{L\omega}\right]=\boldsymbol{\Psi}_\text{Coul}\left(\omega t,f,\frac{g}{L\omega^2},\frac{D}{L},\frac{\rho L^3}{m},\mu,\mu_w\right)
\left[\frac{P}{Mg\sqrt{Lg}},\frac{V}{\sqrt{Lg}}\right]=\boldsymbol{\Psi}_\text{Coul}\left(\sqrt{\frac{g}{L}} t,f,\frac{g}{L\omega^2},\frac{D}{L},\frac{\rho L^3}{M},\mu,\mu_w\right)
\label{CoulscalND}
\end{equation}
%If the gravity, wheel roughness, and granular media are fixed, the values of $g$, $\rho$,  $\mu$, and $\mu_w$ can be absorbed into the function giving the reduced form:
%\begin{equation}
%\left[\frac{P}{mL^2\omega^3},\frac{V}{L\omega}\right]=\boldsymbol{\Psi}_\text{Coul}\left(\omega t,f,\frac{1}{L\omega^2},\frac{D}{L},\frac{L^3}{m}\right)
%\left[\frac{P}{M\sqrt{L}},\frac{V}{\sqrt{L}}\right]=\boldsymbol{\bar{\Psi}}_\text{Coul}\left(\sqrt{\frac{1}{L}} t,f,\frac{1}{L\omega^2},\frac{D}{L},\frac{L^3}{M}\right)
%\label{CoulscalD}
%\end{equation}

This can be further reduced, if it is assumed that granular motion under the wheel is approximately invariant in the out-of-plane dimension.  
In this case, if the mass $M$ and width $D$ of the wheel are scaled by some $C_0$, this would be identical to running $C_0$ copies of the wheel side by side.   
The resulting power would be $C_0 P$ and the velocity would remain unchanged.  From Eq \ref{CoulscalND}, this means $\boldsymbol{{\Psi}}_\text{Coul}$ is unchanged under such a transformation, which 
constrains $\boldsymbol{{\Psi}}_\text{Coul}$ to depend on $M$ and $D$ only through the ratio $D/M$, requiring $\boldsymbol{{\Psi}}_\text{Coul}$ to depend on $D/L$ and $\rho L^3/M$ only through their product:
\begin{equation}
%\left[\frac{P}{mL^2\omega^3},\frac{V}{L\omega}\right]=\boldsymbol{\Psi}_\text{Coul}\left(\omega t,f,\frac{1}{L\omega^2},\frac{D L^2}{m}\right)
\left[\frac{P}{Mg\sqrt{Lg}},\frac{V}{\sqrt{Lg}}\right]=\bar{\boldsymbol{\Psi}}_\text{Coul}\left(\sqrt{\frac{g}{L}} t,f,\frac{g}{L\omega^2},\frac{\rho D L^2}{M},\mu,\mu_w\right)
\label{RFTCoulscalD}
\end{equation}
This form is identical to Eq \ref{RFTscalD} when $g$, $\rho$,  $\mu$, and $\mu_w$ are fixed so that they can be absorbed into the function:
\begin{equation}
%\left[\frac{P}{mL^2\omega^3},\frac{V}{L\omega}\right]=\boldsymbol{\Psi}_\text{Coul}\left(\omega t,f,\frac{1}{L\omega^2},\frac{D L^2}{m}\right)
\left[\frac{P}{M\sqrt{L}},\frac{V}{\sqrt{L}}\right]=\tilde{\boldsymbol{\Psi}}_\text{Coul}\left(\sqrt{\frac{1}{L}} t,f,\frac{1}{L\omega^2},\frac{D L^2}{M}\right)
\label{RFTCoulscalDabs}
\end{equation}
Therefore, the scalings implied by RFT also arise if one assumes the grains obey Coulomb plasticity. This connection presumes the flow under the wheel is invariant in the out-of-plane dimension, per a wheel with large enough $D$ relative to sinkage.  \\

\section{Experiments}
%\noindent{{\textit{Experiments:}} 
To test the above scalings, an experimental sand bed in the MIT Robotic Mobility Group Lab \cite{senatore2014analysis} was used.  The test apparatus, as shown in Fig \ref{ExpSetUp}, consists of a Lexan bin filled with Quikrete medium sand surrounded by an aluminum frame.  Attached to the aluminum frame are two low friction rods, which are attached to the carriage and allow for horizontal motion.  The carriage is also attached through low-friction vertical rods to the main platform to allow for vertical wheel motion.  This allows the wheel to translate freely, e.g. non-circular wheels ``bobble'' in the vertical direction as they drive through the grains.
\begin{figure}[t]
\begin{center}
\includegraphics[width=\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{Fig2_v2.pdf}
%\includegraphics[width=.85\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{Fig2_v2.pdf}
\end{center}
%\includegraphics[width=0.3\textwidth,trim={4.3cm 3.4cm 1.5cm 2.0cm},clip]{LabSetUpAndFBD.pdf}
\caption{Experimental apparatus.  The wheel rotates in the counter-clockwise direction.  A constant force ``spring'' and a pulley are used to alter the wheel's effective gravity.}
%(a)\hspace{-.2cm} \includegraphics[width=0.27\textwidth,trim={4.3cm 3.4cm 1.7cm 2.0cm},clip]{CompiledImages/LabSetUp.pdf} \hspace{-.3cm} (b)\ \ \includegraphics[width=0.16\textwidth,trim={5.8cm 2.4cm 5.7cm 3.1cm},clip]{CompiledImages/FBDcompile.pdf}
%\caption{(a) Experimental apparatus.  The wheel rotates in the counter-clockwise direction.  A constant force ``spring'' and a pulley are used to alter the wheel's effective gravity. (b) Forces on the wheel system. \kcom{I can't see Fres!  Is that easy for color-normal people to see?} \jcom{Does this look better?}}
%\kcom{Should we indicate somewhere, maybe in the text or in one of these images, which way the lug spins? Also, maybe the Fres vector should be slanted %so it's clear it doesn't have to point straight up.}\jcom{Ok added the direction into the caption and fixed the arrow direction.}\kcom{It's a bit weird to show %$F_p$ and $M_p g$ on the same FBD because those would never occur together on one FBD; one belongs to the pulley FBD and the other the the wheel %platform FBD.  Also, we should swap $g$ with $g_{e}$ to be consistent with the below since $g_e$ is the earth value.}}\jcom{Agreed, which one do we %keep though?  If its truly the FBD of the platform it should be $F_p$, but I think the $M_p g_e$ gives more information.  I kept that one but can easily %change.}
\label{ExpSetUp}
\end{figure}
The main platform is connected to a motor driven wheel.  The wheels were printed in PLA plastic using a MakerBot 3D printer.  Position, torque, velocity, and rotational velocity sensors are attached to the frame to record the relevant data.

%Additionally, we developed a setup, as shown in Fig \ref{ExpSetUp}, that allows us to vary the effective gravity on the wheel.  This lets us check Eq \ref{RFTscalND} under different selections of $g$, but keeping $\xi$ held fixed. Lowering gravity also prevents overload of the wheel motor as mass increases.   An SDP/SI Neg'ator constant-force ``spring'' ($F_c=66.7$N) is attached between the carriage and the main platform.  A pulley is attached to the carriage platform and a belt around the pulley is attached to the main platform and a mass, $M_p$.  
Additionally, we developed a setup (Fig \ref{ExpSetUp}) that allows us to vary the effective gravity on the wheel.  This lets us check Eq \ref{RFTscalND} under different selections of $g$, while keeping $\xi$ held fixed. Lowering gravity on the wheel also prevents overload of the wheel motor as mass increases.   An SDP/SI Neg'ator constant-force ``spring'' ($F_c=66.7$N) is attached between the carriage and the main platform.  A pulley is attached to the carriage platform and a belt around the pulley is attached to the main platform and a mass, $M_p$.  In view of Fig 2,  the gravity change can be seen by applying Newton's second law to the main carriage coupled to the hanging weight.  The result is that the wheel's translational motion always matches that of a free wheel whose mass is $M=M_T+M_p$ and gravity is:
\begin{equation}
g=\frac{M_Tg_e-M_pg_e-F_c}{M_T+M_p}
\label{PulleyMass}
\end{equation}
where $g_e$ is earth gravity and $M_T$ is the total mass of the wheel, motor, main platform, and added mass.  This conversion ensures the wheel's gravitational and inertial masses always agree but does not change the gravity experienced by the grains; the grains are still loaded by earth gravity. Hence, the $\xi$ value remains fixed. Consequently, when we use the conversion technique above to relate two tests, the same effective $g$ must be used for both tests. The upcoming section on Potential Extraplanetary Applications uses discrete element simulations, which allow us to test scalings that arise as the true ambient gravity, acting both on the grains and the wheel, is varied between pairs of tests. 
%\begin{figure}[t!!!!!!!!!!!!!!!!!!!]
%\includegraphics[width=0.2\textwidth,trim={5.8cm 2.9cm 5.7cm 3.1cm},clip]{CompiledImages/FBDcompile.pdf}
%\caption{Free body diagram of the spring-pulley system, where the pulley and constant force "spring" are assumed frictionless.  The pulley mass $M_p$ and the tire mass $M_T$ are changed to alter the gravity acting on the wheel.}
%\label{ExpSetUpFBD}
%\end{figure}
%The spring and pulley are used in tandem to vary the effective gravity on the wheel in a way that ensures the wheel's gravitational and inertial masses always agree. This does not change the gravity experienced by the grains; hence the $\xi$ value remains fixed.
%The free body diagram of the system is shown in Figure \ref{ExpSetUpFBD}, where the pulley is assumed to be frictionless such that $M_p g=F_p$.  
%In view of Fig \ref{ExpSetUp},  the gravity change can be seen by applying Newton's second law to the main carriage coupled to the hanging weight.  The result is that the wheel's translational motion always matches that of a free wheel whose mass is $M=M_T+M_p$ and gravity is $g=(M_Tg_e-M_pg_e-F_c)/(M_T+M_p)$ where $g_e$ is earth gravity and $M_T$ is the total mass of the wheel, motor, main platform, and added mass. 
%\begin{equation}
%\begin{bmatrix}F_{\text{res}_x} \\ F_{\text{res}_y}\end{bmatrix} + \frac{M_T-M_p-\frac{F_c}{g_{e}}}{K}\begin{bmatrix}0 \\ -Kg_{e}\end{bmatrix}=\left(M_T+M_p\right)\begin{bmatrix}\ddot x \\ \ddot y\end{bmatrix}
%\label{SumFyKg}
%\end{equation}
%where ${F}_{\text{res}}$ is the resistive force, $F_c$ is the constant force from the spring, $M_p$ is the mass hung from the pulley, and $M_T$ is the mass of the wheel-plus-platform.  The above system is identical to one in which the inertial and gravitational masses are both identically $M=M_T+M_p$ and $g=Kg_e=\frac{M_Tg_e-M_pg_e-F_c}{M_T+M_p}$.
%\begin{equation}
%\frac{M_T-M_p-\frac{F_c}{g}}{K}=\left(M_T+M_p\right)
%\label{GravInert}
%\end{equation}

The scalings from Equation \ref{RFTscalND} were tested with a set of 288 experiments involving pairs of lug or cylindrical wheels (Fig \ref{WheelPairs}) having different size dimensions, mass loadings, and rotation speeds \footnote{The interior circle on the lug wheels is for mounting and never came into contact with the sand in any tests.}.  \edit{Cylindrical wheels are always} covered in sandpaper to increase interface friction.  Lug wheels have four arms and an elbow half way down their lengths bent 150$^{\circ}$\edit{; previous work suggests an elbow bend improves wheel efficiency \cite{li13}.  The two wheel shapes (cylindrical and lug) were chosen to demonstrate the scaling over two distinct driving motions; cylindrical wheels drive smoothly with limited sinkage, while the lug wheels dig through and remove pockets of sand as they translate (see supplemental movies).} 
%\begin{figure}
%\includegraphics[width=0.6\columnwidth,trim={0cm 0cm 0cm 0cm},clip]{Cyl4Bar.pdf}
%\caption{The two shapes, $f$, used in our study: cylinders (left) and lugs (right), and their corresponding definitions of $L$.  The interior circle on the lug wheels was used for mounting and never came into contact with the sand.}
%\label{Cyl4bar}
%\end{figure}

\begin{table}
\caption{Experimental tests performed and inputs used.  Each row represents a different wheel pair for which 12 combinations of driving parameters were applied; mass and rotational velocity pairs represented as cartesian products.}
\begin{center}
\begin{tabular}{| >{\centering\arraybackslash} m{0.6cm}| >{\centering\arraybackslash} m{3.8cm}| >{\centering\arraybackslash} m{3.7cm}|}
\hline
 & Common Parameters \: [cm], [cm], [m/s$^2$], [ ] & Driving Parameters \newline [kg] x [deg/s] \\
\hline
Key & $\scriptscriptstyle\begin{bmatrix} L_b\\ L_s \end{bmatrix}$,$\scriptscriptstyle\begin{bmatrix} D_b\\ D_s \end{bmatrix}$,$\scriptscriptstyle\begin{bmatrix} g \end{bmatrix}$,$\displaystyle\left(\{\text{shape}\}\right)$  & \:\: {$\scriptscriptstyle\left(\begin{bmatrix} M_b^1\\ M_s^1 \end{bmatrix}\begin{bmatrix} M_b^2\\ M_s^2 \end{bmatrix}\begin{bmatrix} M_b^3\\ M_s^3 \end{bmatrix}\right)$ x \:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\: $\scriptscriptstyle\left(\begin{bmatrix} \omega_b^1\\ \omega_s^1 \end{bmatrix}\begin{bmatrix} \omega_b^2\\ \omega_s^2 \end{bmatrix}\begin{bmatrix} \omega_b^3\\ \omega_s^3 \end{bmatrix}\begin{bmatrix} \omega_b^4\\ \omega_s^4 \end{bmatrix}\right)$} \:\:\:\:\:\:\:\:\:\:\:\: \\
\hline
Pair A & $\scriptscriptstyle\begin{bmatrix} 12.50\\ 8.08 \end{bmatrix}$,$\scriptscriptstyle\begin{bmatrix} 15\\ 15 \end{bmatrix}$,$\scriptscriptstyle\begin{bmatrix} 3.71 \end{bmatrix}$,$\displaystyle\left(\{\text{Cyl}\}\right)$  & \:\:{$\scriptscriptstyle \left(\begin{bmatrix} 35.9\\ 14.9 \end{bmatrix}\begin{bmatrix} 42.2\\ 17.6 \end{bmatrix}\begin{bmatrix} 45.7\\ 19.0 \end{bmatrix}\right)$ x \:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\: $\scriptscriptstyle\left(\begin{bmatrix} 14.0\\ 17.4 \end{bmatrix}\begin{bmatrix} 17.0\\ 21.2 \end{bmatrix}\begin{bmatrix} 20.0\\ 24.9 \end{bmatrix}\begin{bmatrix} 23.0\\ 28.6 \end{bmatrix}\right)$} \:\:\\
\hline
Pair B & $\scriptscriptstyle\begin{bmatrix} 11.25\\ 7.50 \end{bmatrix}$,$\scriptscriptstyle\begin{bmatrix} 14\\ 14 \end{bmatrix}$,$\scriptscriptstyle\begin{bmatrix} 1.31 \end{bmatrix}$,$\displaystyle\left(\{\text{Lug}\}\right)$  & \:\:{$\scriptscriptstyle\left(\begin{bmatrix} 30.2\\ 13.4 \end{bmatrix}\begin{bmatrix} 34.9\\ 15.5 \end{bmatrix}\begin{bmatrix} 39.7\\ 17.6 \end{bmatrix}\right)$ x \:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\: $\scriptscriptstyle \left(\begin{bmatrix} 14.0\\ 17.1 \end{bmatrix}\begin{bmatrix} 17.0\\ 20.8 \end{bmatrix}\begin{bmatrix} 20.0\\ 24.5 \end{bmatrix}\begin{bmatrix} 23.0\\ 28.2 \end{bmatrix}\right)$} \:\: \\
\hline
Pair C & $\scriptscriptstyle\begin{bmatrix} 11.25\\ 9.00 \end{bmatrix}$,$\scriptscriptstyle\begin{bmatrix} 14\\ 10 \end{bmatrix}$,$\scriptscriptstyle\begin{bmatrix} 1.31 \end{bmatrix}$,$\displaystyle\left(\{\text{Lug}\}\right)$  & \:\: {$\scriptscriptstyle\left(\begin{bmatrix} 29.3\\ 13.4 \end{bmatrix}\begin{bmatrix} 34.0\\ 15.5 \end{bmatrix}\begin{bmatrix} 38.6\\ 17.6 \end{bmatrix}\right)$ x \:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\: $\scriptscriptstyle\left(\begin{bmatrix} 14.0 \\ 15.6 \end{bmatrix}\begin{bmatrix} 17.0\\ 19.0 \end{bmatrix}\begin{bmatrix} 20.0\\ 22.4 \end{bmatrix}\begin{bmatrix} 23.0\\ 25.7 \end{bmatrix}\right)$} \:\: \\
\hline
\end{tabular}
\end{center}
\label{ExpInputs}
\end{table}

\begin{figure}[t!!!!!!!!!!!!!!!!!!!!]
(a) %\includegraphics[width=0.6\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{PairAWheels.png}
\includegraphics[width=0.6\columnwidth,trim={2.0cm 9.3cm 2.7cm 8.5cm},clip]{PairAWheels.pdf} 
\\
\vspace{0.05cm}
(b)  \includegraphics[width=0.6\columnwidth,trim={2.0cm 8cm 2.0cm 8cm},clip]{PairBWheels.pdf} 
\\
\vspace{0.05cm}
(c) \includegraphics[width=0.6\columnwidth,trim={2.0cm 8.5cm 2.0cm 8.5cm},clip]{PairCWheels.pdf}
\caption{The three different wheel pairs tested: Pair A (a), Pair B (b), and Pair C (c).}
\label{WheelPairs}
\end{figure}

The different test inputs are shown in Table \ref{ExpInputs}. Four repetitions of the same test are performed each time to obtain useful error-bars. In each test, a pair of wheels of the same shape but different size dimensions were driven through the sand; one wheel is `big' and one is `small' as denoted with subscripts $b$ and $s$ in the table.  Three pairs of wheels were used in total as shown in Fig \ref{WheelPairs}.  In each test, the two masses ($M_b$ and $M_s$) and rotation speeds ($\omega_b$ and $\omega_s$) were chosen to ensure the big and small systems have the same dimensionless inputs to $\boldsymbol{\Psi}$.  The time-averaged non-dimensional power $\langle \tilde{P} \rangle=\langle \frac{P}{Mg\sqrt{Lg}}\rangle$ and velocity $\langle \tilde{V} \rangle=\langle \frac{V}{\sqrt{Lg}}\rangle$ were measured when the wheel motion reached cyclic behavior.  \edit{The power was calculated by multiplying the constant rotation speed by the average torque measured over one cycle.  The velocity was measured directly and averaged over one cycle.  The wheels were started in the same position, such that the torque and velocity measurements could be averaged over corresponding cycles.}  In each case,  the proposed scaling relation is satisfied if $\langle \tilde{P}_b \rangle=\langle \tilde{P}_s \rangle$ and $\langle \tilde{V}_b \rangle=\langle \tilde{V}_s \rangle$. 

The first tests performed were with cylindrical wheels (Pair A).  The non-dimensional powers and velocities arising from each of the 12 pairs of driving parameters are plotted in Fig \ref{alldata}(a).  In general the expected trend is observed at low rotational velocity, with some deviation at the highest rotational rate.
\begin{figure}[t!!!!!!!!!!!!!!!!!!!!]
\includegraphics[width=\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{fig_combined.pdf}
%(a) \includegraphics[width=0.75\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{CylTrial.pdf} 
%\\
%\vspace{0.05cm}
%(b)  \includegraphics[width=0.75\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{SlapTrial.pdf} 
%%(b) \includegraphics[width=0.9\columnwidth,trim={2.1cm 1.2cm 1.2cm 0.5cm},clip]{Pictures/SlapNDPV2.pdf} 
%\\
%\vspace{0.05cm}
%(c) \includegraphics[width=0.75\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{All4Trial.pdf}


%(a) \includegraphics[width=0.75\columnwidth,angle=90,trim={0.0cm 0.7cm 0.0cm 0.5cm},clip]{CylNDPV3.pdf} 
%\\
%\vspace{0.05cm}
%(b) \includegraphics[width=0.75\columnwidth,angle=90,trim={0.0cm 0.7cm 0.0cm 0.9cm},clip]{SlapNDPV4.pdf} 
%%(b) \includegraphics[width=0.9\columnwidth,trim={2.1cm 1.2cm 1.2cm 0.5cm},clip]{Pictures/SlapNDPV2.pdf} 
%\\
%\vspace{0.05cm}
%(c) \includegraphics[width=0.75\columnwidth,angle=90,trim={0.0cm 0.8cm 0.0cm 0.8cm},clip]{All4NDPV3.pdf}
\caption{Time-average non-dimensional power and velocity values for {Pair A} (a), {Pair B} (b), and {Pair C} (c).  Each color in each plot represents a different mass pair selection and the corresponding four points of each color represent different rotational velocity pairings.  The solid line is the model prediction of $\frac{\langle\tilde{P}_b\rangle}{\langle\tilde{P}_s\rangle}=\frac{\langle\tilde{V}_b\rangle}{\langle\tilde{V}_s\rangle}=1$.
%In this case, displaying that the non-dimensional power and velocity of each wheel are close to equivalent represents a scaling of the \textbf{dimensional} power and velocity as $P_b=2.99P_s$ and $V_b=1.25V_s$.
}
\label{alldata}
%\label{CylNDPV}
\end{figure}
The next series of tests utilized two lug wheels (Pair B).  %The average non-dimensional power and velocity over a single cycle are plotted in Figure \ref{SlapNDPV}.  
Again, the expected scaling relation is observed (Fig \ref{alldata}(b)).
%\begin{figure}[h!!!!!!!!!!!!!!!!!!!!]
%\includegraphics[width=0.99\columnwidth,trim={2.1cm 1.2cm 1.2cm 0.5cm},clip]{Pictures/SlapNDPV2.pdf} 
%\caption{\textbf{Pair B} - Average non-dimensional power and velocity values over a single cycle for each lug wheel pair are plotted.  Each color represents a different mass pair and the corresponding four points of each color represent a different rotational velocity pair. 
%% In this case, displaying that the non-dimensional power and velocity of each wheel are close to equivalent represents a scaling of the %\textbf{dimensional} power and velocity as $P_b=2.76P_s$ and $V_b=1.22V_s$.
% }
%\label{SlapNDPV}
%\end{figure}
The final series of experiments performed (Pair C) also used four-arm lug wheels, however the widths and lengths of the two wheels were not equal and not proportional.  Hence, each test in this series involves two wheels with different lengths, widths, masses, and rotation speeds --- all differing by different factors --- making this series the most stringent, and arguably most interesting, test of the proposed scaling relation. Strong agreement is observed (Fig \ref{alldata}(c)). The best-fit slopes of all six datasets in Fig \ref{alldata} are all within $3\%$ of the predicted value of 1. %\edit{Further information is provided in the Supplemental Materials regarding how the actual time-dependence of these quantities can also be related.}
%of the form $\langle \tilde{P}_b \rangle=A\langle \tilde{P}_s \rangle$ and $\langle \tilde{V}_b \rangle=A\langle \tilde{V}_s \rangle$ were found to be within 3\% of the predicted slope of 1.
%\begin{figure}[h!!!!!!!!!!!!!!!!!!!!]
%\includegraphics[width=0.99\columnwidth,trim={2.0cm 1.2cm 1.2cm 0.4cm},clip]{Pictures/All4NDPV2.pdf} 
%\caption{\textbf{Pair C} - Average non-dimensional power and velocity values over a single cycle for each lug wheel pair are plotted.  Each color represents a different mass pair and the corresponding four points of each color represent a different rotational velocity pair.  
%%In this case, displaying that the non-dimensional power and velocity of each wheel are close to equivalent represents a scaling of the 
%%\textbf{dimensional} power and velocity as $P_b=2.45P_s$ and $V_b=1.12V_s$.
%}
%\label{All4NDPV}
%\end{figure}
% The average non-dimensional power and velocity over a single cycle for all mass and rotational velocity pairs are plotted in Figure \ref{All4NDPV}. 


%From Eq \ref{RFTscalND}, in addition to the time-average behaviors, the actual time-dependence of these quantities should relate when time is accordingly scaled.  Figure \ref{SlapScalePowUp} shows the non-dimensional power plotted against dimensionless time, $\tilde t = \sqrt{\frac{g}{L}}t$, over the course of a single cycle from two pairs of experiments from wheel Pair C.  The predicted \edit{scaling} agreement between trajectories is observed.  The result is non-trivial; the dimensional powers differ by a factor of 2.45 for each pair. 

\edit{In mitigating several laboratory constraints on the size, speed, and mass loading of our wheels, we selected the range of input variables to ensure an ample variation of all inputs rather than a large range of any one variable.   This was done to ensure the dimensionless groups we have identified robustly describe the locomotion process, and that any additional dimensionless groups, which more complex sand models might incur, are not crucial to describe the driving outputs.}

%\section{Relating Time-Dependent Parameters} 
From Eq 2, in addition to the time-average behaviors, the actual time-dependence of $\tilde P$ and $\tilde V$ should relate when time is accordingly scaled.  Figure \ref{SlapScalePowUp} shows the non-dimensional power plotted against dimensionless time, $\tilde t = \sqrt{\frac{g}{L}}t$, over the course of a single cycle from two pairs of experiments from wheel Pair C.  The predicted \edit{scaling} agreement between trajectories is observed.  The result is non-trivial; the dimensional powers differ by a factor of 2.45 for each pair. 

\begin{figure}[h!!!!!!!!!!!!!!!!!!!!!!!!!!!!]
%\includegraphics[width=\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{Fig3.png}
\includegraphics[width=0.55\columnwidth,angle=-90,trim={4.0cm 4.0cm 4.0cm 4.0cm},clip]{Fig3CombSS.pdf}
%\includegraphics[width=0.96\columnwidth,trim={2.1cm 1.6cm 3.3cm 2.3cm},clip]{UpdatedNDPower5.pdf}
%\includegraphics[width=0.75\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{UpdatedNDPowerTrial.pdf}
\caption{Non-dimensional power trajectories for two pairs of lug-wheel experiments over a single cycle \edit{(90$^{\circ}$ of rotation)}.  One pair is shown in solid lines, while the other is shown in dotted lines.  }
\label{SlapScalePowUp}
\end{figure}
%  For additional generality, recall that our scaling law derives from both Coulomb-Plasticity and RFT, hence it should work if either theory is applicable.  Previous research suggests these theories accurately describe a broad range of problems ---  Coulomb plasticity has been shown to model slow and fast flows in diverse geometries \cite{dunatunga2015continuum,askari2015intrusion,lagree2011granular} and RFT has shown valid in much faster locomotion problems involving robots running many times their body-length per second  \cite{li13}.}\\
%\begin{figure}[h!!!!!!!!!!!!!!!!!!!!!!!!!!!!]
%\includegraphics[width=0.85\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{Fig3Comb.pdf}
%\includegraphics[width=0.96\columnwidth,trim={2.1cm 1.6cm 3.3cm 2.3cm},clip]{UpdatedNDPower5.pdf}
%\includegraphics[width=0.75\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{UpdatedNDPowerTrial.pdf}
%\caption{Non-dimensional power trajectories for two pairs of lug-wheel experiments over a single cycle \edit{(90$^{\circ}$ of rotation)}.  One pair is shown in solid lines, while the other is shown in dotted lines.  }
%\label{SlapScalePowUp}
%\end{figure}

\section{Potential Extraplanetary Applications}
%\noindent{\textit{Potential extraplanetary applications:}} 
%Our experimental apparatus allows the gravity acting on the wheel to be reduced, but leaves $\xi$ unchanged.  This did not hinder our ability to validate Eq \ref{RFTscalND} because each test pair considered had common $\xi$ and $g$ values.   However,
%It would be advantageous to be able to scale to two experiments in the same sand but with different ambient gravities; e.g. to use an earthbound experiment to predict an extraplanetary experiment. This would require a theory to explain how $\xi$ varies with gravity in Eq \ref{RFTscalND}.  
It would be advantageous to be able to scale to two experiments in the same sand but with different ambient gravities; e.g. to use an earthbound experiment to predict an extraplanetary experiment. To derive the proper scalings, we can use Eq \ref{RFTCoulscalD}, for which gravity dependence of the granular response is natively accounted.  Interestingly, we can show that the same scalings arise assuming RFT (Eq \ref{RFTscalND}) once the dependence of $\xi$ on gravity is deduced.  %\edit{Further analysis is provided in the Supplemental Materials to describe scaling relations when the wheel is also horizontally loaded with a backward drawbar force.} \\
To deduce a form for $\xi$, we assume the grains behave as an ideal Coulomb material, and then dimensional analysis implies the formula for $\xi$ obeys
\begin{equation}
\xi=\rho g \hat\xi\left(\mu,\mu_w\right).
\label{xig}
\end{equation}
  %Alternatively, we can simply use Eq \ref{RFTCoulscalD}, for which gravity dependence of the granular response is natively accounted. 

%\begin{figure}[t]
%\includegraphics[width=\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{DEM.pdf}
%\includegraphics[width=0.96\columnwidth,trim={2.1cm 1.6cm 3.3cm 2.3cm},clip]{UpdatedNDPower5.pdf}
%\includegraphics[width=0.75\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{UpdatedNDPowerTrial.pdf}
%\caption{\edit{DEM simulations: Velocity data for three test pairs (six simulations).   The `big' simulation of each pair has $4\times$ the gravity of the `small'.  Selection of $\tilde{\rho}$ and $\tilde{g}$ produces different modes of locomotion (e.g. walking, trudging, skipping).}}\label{DEM}
%\end{figure}

By substituting Eq \ref{xig} into Eq \ref{RFTscalND}, we obtain
\begin{equation}
\left[\frac{P}{Mg\sqrt{Lg}},\frac{V}{\sqrt{Lg}}\right]=\boldsymbol{\Psi}\left(\sqrt{\frac{g}{L}} t,f,\frac{g}{L\omega^2},\frac{\rho \hat\xi\left(\mu,\mu_w\right) D L^2}{M}\right).
%\left[\frac{P}{Mg\sqrt{Lg}},\frac{V}{\sqrt{Lg}}\right]=\boldsymbol{\Psi}\left(\sqrt{\frac{g}{L}}  t,f,\frac{g}{L\omega^2},\frac{\rho \hat\xi\left(\mu,\mu_w\right) D L^2}{M}\right)
\label{RFTgrav}
\end{equation}
Note that this relation is equivalent to Eq \ref{RFTCoulscalD} as long as material parameters for the grains are absorbed into the function.  So, even when gravity varies between tests, both approaches again yield the same scaling relations, which can be expressed as follows: Consider two experiments with common wheel shape $f$ and common grains, where one test is described by the inputs $(g, L, M, D,\omega)$  and the other by the inputs $(g', L',M',D',\omega')=(qg, r L, s M, sr^{-2}D, q^{1/2}r^{-1/2}\omega)$ for any positive scalars $q$, $r$, and $s$. Then the steady driving cycles of the corresponding outputs should obey $\langle P'\rangle=q^{3/2} r^{1/2}s \langle P\rangle$ and $\langle V'\rangle=q^{1/2}r^{1/2}\langle V\rangle$. 


%\edit{
%\noindent{\textit{Towing a drawbar force:}} 
%One useful consideration in wheel design is the ability of a wheel to pull a load acting in the opposite direction of motion.  This can be described as a constant drawbar force $F_d$, which acts in the negative horizontal direction.  With this new consideration, we can extend the relationship found in Eq \ref{RFTCoulscalD} to include $F_d$ by adding an additional non-dimensional group:
%\begin{multline}
%\left[\frac{P}{Mg\sqrt{Lg}},\frac{V}{\sqrt{Lg}}\right]=\\
%\boldsymbol{\Psi}\left(\sqrt{\frac{g}{L}}  t,f,\frac{g}{L\omega^2},\frac{\rho \hat\xi\left(\mu,\mu_w\right) D L^2}{M},\frac{F_d}{Mg}\right)
%\label{RFTCoulscalD}
%\end{multline}
%Using this new form, we can expand the scaling law to include the drawbar force.  Again, both experiments would have common $f$, $\rho$, $\mu$, and $\mu_w$, but one is described by the inputs $(g, L, M, D,\omega, F_d)$  and the other by the inputs $(g', L',M',D',\omega', F_d')=(qg, r L, s M, sr^{-2}D, q^{1/2}r^{-1/2}\omega, sq F_d)$ for any positive scalars $q$, $r$, and $s$. The steady driving cycles of the corresponding outputs should then again obey $\langle P'\rangle=q^{3/2} r^{1/2}s \langle P\rangle$ and $\langle V'\rangle=q^{1/2}r^{1/2}\langle V\rangle$.
%}\\

\begin{figure}
\begin{center}
(a)\ \ \ \ \ \ \includegraphics[width=.5\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{walk_pic.pdf}
\\
(b)\includegraphics[width=.6\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{walk_velocity.pdf}
\\
(c)\includegraphics[width=.6\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{walk_power.pdf}
\end{center}
\caption{\textbf{Walking} mode ($\tilde{\rho}=44.4$, $\tilde{g}=14.8$).  Small test inputs: $g_s=9.8\frac{\text{m}}{\text{s}^2}$, $L_s= .0168\text{m}$, $M_s=12.7\times10^{-6}\text{kg}$, $\omega_s=6.28 \frac{\text{rad}}{\text{s}}$.  Big test inputs: $g_b=39.2\frac{\text{m}}{\text{s}^2}$, $L_b= .0224\text{m}$, $M_b=22.4\times10^{-6}\text{kg}$, $\omega_b=10.9 \frac{\text{rad}}{\text{s}}$. Both tests use the same grains.  (a) Snapshot of the small and big systems at a common value of $\tilde{t}$. (b)-(c)  Time-dependence of power and velocity in dimensionless form during locomotion process. Blue lines are for the big test and red lines for the small.}\label{walk}
\end{figure}

\edit{To see if indeed the above form can be used to relate to tests in two different gravities, we perform a sequence of simulations using Discrete Element Method (DEM). DEM tests allow for the easy manipulation of input parameters, including gravitational acceleration $g$, over a broad scope of values not limited by experimental hardware. We use LAMMPS \cite{plimpton1995fast} for the simulations. Throughout, the grains have the following characteristics: density $\rho_s = 2500$ $\frac{\text{kg}}{\text{m}^3}$, normal and tangential stiffnesses $k_N = 16411 \frac{\text{N}}{\text{m}}$ and $k_t = 4689 \frac{\text{N}}{\text{m}}$, normal damping coefficients $\gamma_N = 184835 \text{s}^{-1}$ and friction coefficient $\mu = 0.4$. The simulations are effectively 2D; the bed is modeled as a planar packing of polydisperse spheres with a linearly distributed radius with minimum  $R_{\text{min}} = 0.254$mm and maximum  $R_{\text{max}} = 0.381$mm. These choices for the particle material have been used successfully in past studies \cite{koval2009annular} and resemble the parameters for glass beads. The wheels are made of a coplanar collection of $0.8$mm diameter particles constrained to move in the plane as a rigid body ($D=0.8$mm).

%\begin{tabular}{|l|l|}0
%\hline
%  &  Wheel parameters \:  [m/s$^2$], [m], [kg], [rad/s ] 
 % \\
%\hline
%Key &  $\scriptscriptstyle\begin{bmatrix} g_b \\ g_s \end{bmatrix}$,\ $\scriptscriptstyle\begin{bmatrix} L_b\\ L_s \end{bmatrix}$,\  $\scriptscriptstyle\begin{bmatrix} M_b\\ M_s \end{bmatrix}$,\ $\scriptscriptstyle\begin{bmatrix} \omega_b\\ \omega_s \end{bmatrix}$\\
%\hline
%Walk &  $\scriptscriptstyle\begin{bmatrix} 39.2 \\ 9.8 \end{bmatrix}$,\ $\scriptscriptstyle\begin{bmatrix}0.0224 \\ 0.0168 \end{bmatrix}$,\  $\scriptscriptstyle\begin{bmatrix} 2.24e-5\\ 1.27e-5 \end{bmatrix}$,\ $\scriptscriptstyle\begin{bmatrix} 10.9 \\ 6.28 \end{bmatrix}$\\
%\hline
%end{tabular}


Per Eq \ref{RFTCoulscalD}, when material is fixed, the key inputs are  $\tilde{\rho}\equiv\frac{\rho D L^2}{M}$, $\tilde{g}  \equiv \frac{g}{L\omega^2}$, and wheel shape $f$.  All DEM tests performed use $f$ corresponding to a fully rough bar rotating about its center, with $L$ the length of the bar.  The scaling law is checked by running pairs of simulations with common values of $\tilde{\rho}$ and $\tilde{g}$. In each pair, the dimensional inputs for the `big' test are given from the `small' by $(g_b, L_b, M_b, \omega_b)=(4g_s,\frac{4}{3}L_s,\frac{16}{9}M_s,\sqrt{3}\omega_s)$. Choosing the factors this way guarantees each of the pair has the same dimensionless inputs according to Eq \ref{RFTCoulscalD}. Since all four inputs differ in each test pair, this offers a stringent test of the scaling relation. We emphasize that the big wheel test has $4\times$ the ambient gravity of the small in each test pair, and the same grains are used among all tests.
%for each wheel, but having the same dimensionless parameters $\tilde{\rho}$ and $\tilde{g}$. However, the members of each pair contain very different dimensional parameters, with one 'small' wheel which experiences earth gravity and one 'big' wheel which experiences four times as much gravity. The dimensional parameters that are varied each test are $(g, L, M, \omega)$. 

The first pairing has small wheel inputs as follows: gravity is $g_s=9.8\frac{\text{m}}{\text{s}^2}$, wheel length is $L_s= .0168\text{m}$, wheel mass is $M_s=12.7\times10^{-6}\text{kg}$, and rotation speed is $\omega_s=6.28 \frac{\text{rad}}{\text{s}})$. The big wheel inputs are $4g_s, \frac{4}{3}L_s,\frac{16}{9}M_s,\sqrt{3}\omega_s$ as noted above. These tests both have $\tilde{\rho}=44.4$, $\tilde{g}=14.8$, and produce a `walking' motion where the wheels are so light compared to the grains they barely deform the granular surface (see supplemental videos).  Figure \ref{walk} displays the results. The scaling law would predict the two tests should have the same $\tilde{V}=V/\sqrt{Lg}$ vs $\tilde{t}$ and $\tilde{P}=P/Mg\sqrt{Lg}$ vs $\tilde{t}$ curves because they share all the same dimensionless inputs to Eq \ref{RFTCoulscalD}. The strong overlap of the $\tilde{V}=V/\sqrt{Lg}$ vs $\tilde{t}=t\sqrt{g/L}$ curves verifies the scaling relation in its full time-dependent form. The snapshot in Fig \ref{walk}(a) shows the two wheels at a common value of $\tilde{t}$. The similarity in the wheel positioning and granular free surface between properly scaled tests at the same $\tilde{t}$ is clear.




The second pairing increases both wheel masses ten-fold, giving $\tilde{\rho}=4.44$, $\tilde{g}=14.8$, which produces a `trudging' motion where the wheels are so heavy they submerge nearly to the central axle while driving through the grains (see supplemental videos).  Results of these tests are presented in Fig \ref{trudge}, which confirms the scaling law's prediction that these two tests should have the same dimensionless power behavior and dimensionless velocity behavior.  Moreover, the snapshot in Fig \ref{trudge}(a) shows that the divit in the free surface created by each wheel looks geometrically similar to the other, which supports the notion that these two problems obey a proper scaling.


\begin{figure}
\begin{center}
(a)\ \ \ \ \ \ \includegraphics[width=.5\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{trudge_pic.pdf}
\\
(b)\includegraphics[width=.6\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{trudge_velocity.pdf}
\\
(c)\includegraphics[width=.6\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{trudge_power.pdf}
\end{center}
\caption{\textbf{Trudging} mode ($\tilde{\rho}=4.44$, $\tilde{g}=14.8$). Small test inputs: $g_s=9.8\frac{\text{m}}{\text{s}^2}$, $L_s= .0168\text{m}$, $M_s=127\times10^{-6}\text{kg}$, $\omega_s=6.28 \frac{\text{rad}}{\text{s}}$.  Big test inputs: $g_b=39.2\frac{\text{m}}{\text{s}^2}$, $L_b= .0224\text{m}$, $M_b=224\times10^{-6}\text{kg}$, $\omega_b=10.9 \frac{\text{rad}}{\text{s}}$. Both tests use the same grains. (a) Snapshot of the small and big systems at a common value of $\tilde{t}$. (b)-(c)  Blue lines are for the big test and red lines for the small.}\label{trudge}
\end{figure}


The third simulation pair uses an intermediate pair of masses ($M_s=42.2\times10^{-6}\text{kg}$) but increases the rotation speeds ten-fold, giving $\tilde{\rho}=13.4$, $\tilde{g}=0.148$.  These inputs cause a `skipping' motion where the wheels lift off the bed entirely each cycle and eject grains as they propel forward (see supplemental videos).  Figure \ref{skip} displays the results for the pair of tests.  There is more noise in the data due to the rapid impulses, but the data still indicates a strong agreement between the big and small tests.  Unlike the prior two cases, the skipping case takes a significant time to accelerate to the steady driving motion; acceleration is evident from $\tilde{t}=0$ to $\sim 15$ in Fig \ref{skip}(b).  The snapshot shows that at a common value of $\tilde{t}$ various observations about the two systems look very similar; the free surface shape and ejected grain plume for example look very similar as one might expect for two problems that are properly scaled versions of each other.


\begin{figure}
\begin{center}
(a)\ \ \ \ \ \ \includegraphics[width=.7\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{skip_pic.pdf}
\\
(b)\includegraphics[width=.6\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{skip_velocity.pdf}
\\
(c)\includegraphics[width=.6\columnwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{skip_power.pdf}
\end{center}
\caption{\textbf{Skipping} mode ($\tilde{\rho}=13.4$, $\tilde{g}=0.148$). Small test inputs: $g_s=9.8\frac{\text{m}}{\text{s}^2}$, $L_s= .0168\text{m}$, $M_s=42.2\times10^{-6}\text{kg}$, $\omega_s=62.8 \frac{\text{rad}}{\text{s}}$.  Big test inputs: $g_b=39.2\frac{\text{m}}{\text{s}^2}$, $L_b= .0224\text{m}$, $M_b=75.0\times10^{-6}\text{kg}$, $\omega_b=109 \frac{\text{rad}}{\text{s}}$. Both tests use the same grains.  (a) Snapshot of the small and big systems at a common value of $\tilde{t}$. (b)-(c)  Blue lines are for the big test and red lines for the small. }\label{skip}
\end{figure}

The tests we have conducted cover a fairly wide range of qualitative behaviors of the locomotor, and it is reassuring that even as the mode of locomotion changes, the scaling relation still appears to hold, as is evident in Figs \ref{walk}-\ref{skip}.
%except in the second pairing each wheel was ten times heavier and in the third pairing each wheel rotated ten times faster.
}

%\section{Towing a Drawbar Force}
%\noindent{\textit{Towing a drawbar force:}} 
%Another useful consideration in wheel design is the ability of a wheel to pull a load which causes force in the opposite direction of motion. This can be described as a constant drawbar force $F_d$, which acts in the backward (negative) horizontal direction.  The drawbar force is useful as an indicator of both tractive performance and ability to scale inclines. With this new consideration, we can extend the relationship found in Eq \ref{RFTgrav} to include $F_d$ by adding an additional non-dimensional group:
%\begin{multline}
%\left[\frac{P}{Mg\sqrt{Lg}},\frac{V}{\sqrt{Lg}}\right]=\\
%\boldsymbol{\Psi}\left(\sqrt{\frac{g}{L}}  t,f,\frac{g}{L\omega^2},\frac{\rho \hat\xi\left(\mu,\mu_w\right) D L^2}{M},\frac{F_d}{Mg}\right)
%\label{RFTCoulscalFD}
%\end{multline}
%Using this new form, we can expound the scaling.  Suppose two experiments with common $f$, $\rho$, $\mu$, and $\mu_w$, but one is described by the inputs $(g, L, M, D,\omega, F_d)$  and the other by the inputs $(g', L',M',D',\omega', F_d')=(qg, r L, s M, sr^{-2}D, q^{1/2}r^{-1/2}\omega, sq F_d)$ for any positive scalars $q$, $r$, and $s$. The steady driving cycles of the corresponding outputs should then obey $\langle P'\rangle=q^{3/2} r^{1/2}s \langle P\rangle$ and $\langle V'\rangle=q^{1/2}r^{1/2}\langle V\rangle$.\\

%\vspace{1em}


\section{Discussion}
%\noindent{\textit{Discussion:}} 
We have proposed and validated a set of scaling laws in granular locomotion for the case of rigid, arbitrarily shaped wheels, which was initially obtained by analyzing RFT system dynamics.  The scaling analysis has been reconciled with Coulomb plasticity and an extension has been \edit{shown} to relate locomotion processes in different ambient gravity.  Our study has shown the scaling relations to work over three different wheel shapes (cylinders, lugs, and bars), and two different choices for the granular material (sand in the experiments, and polydisperse spheres in the DEM tests).

In the future, we would like to extend this scaling to include a wheel's ability to pull a load. This can be described as a constant drawbar force $F_d$, which acts in the backward (negative) horizontal direction.  The drawbar force is useful as an indicator of both tractive performance and ability to scale inclines. With this new consideration, we can extend the relationship found in Eq \ref{RFTgrav} to include $F_d$ by adding an additional non-dimensional group:
\begin{multline}
\left[\frac{P}{Mg\sqrt{Lg}},\frac{V}{\sqrt{Lg}}\right]=\\
\boldsymbol{\Psi}\left(\sqrt{\frac{g}{L}}  t,f,\frac{g}{L\omega^2},\frac{\rho \hat\xi\left(\mu,\mu_w\right) D L^2}{M},\frac{F_d}{Mg}\right)
\label{RFTCoulscalFD}
\end{multline}
Using this new form, we can expound the scaling.  Suppose two experiments with common $f$, $\rho$, $\mu$, and $\mu_w$, but one is described by the inputs $(g, L, M, D,\omega, F_d)$  and the other by the inputs $(g', L',M',D',\omega', F_d')=(qg, r L, s M, sr^{-2}D, q^{1/2}r^{-1/2}\omega, sq F_d)$ for any positive scalars $q$, $r$, and $s$. The steady driving cycles of the corresponding outputs should then obey $\langle P'\rangle=q^{3/2} r^{1/2}s \langle P\rangle$ and $\langle V'\rangle=q^{1/2}r^{1/2}\langle V\rangle$.  It would be interesting future work to experimentally validate this proposed scaling.

In general, there are two possible explanations for why scaling laws arising from rate-independent Coulomb plasticity (or RFT) work as well as they do, even though granular material is known to obey more complicated constitutive behaviors (as described in the Introduction).  First, it is possible that the complexities of granular rheology influence granular kinematics moreso than boundary stresses. Findings in \cite{henann2013predictive} showed that in split-bottom Couette flow, even though the flow features are influenced noticeably by nonlocality due to grain size, the stresses transmitted to the solid boundaries are largely unaffected by nonlocal effects.  This suggests that even if the flow fields predicted by Coulomb plasticity are less accurate than those of more complete models, the intruder forces it predicts may be quite similar.  If so, the intruder dynamics would still be accurately modeled despite using a simplified granular continuum model.  

The second explanation is that locomotion in granular beds, even over a relatively large range inputs, does not bring about the circumstances where additional constitutive phenomena become relevant in the flow.  For example, Coulomb theory neglects rate-sensitivity of the material, which is known to exist (i.e. the $\mu(I)$ rheology \cite{dacruz05,jop06}).  One might expect that for more rapidly spinning wheels this effect could be significant, but oftentimes modifications for rate-dependency have a minimal effect on the flow solution \cite{lagree2011granular} because the rate-effects are scaled by pressure through the inertial number.  %For example,  robots that run on sand using ``c-legs'' \cite{li13} move many times their body-length per second yet remain well-described by rate-independent RFT.  
Size-effects and history-dependent phenomena (e.g. stress-dilatancy and anisotropy evolution) could also be less significant in locomotive intrusion processes since we deal with  macroscopic intruders (as opposed to those competing with the grain size \cite{goldsmith2013drag}) that quickly induce large strains of the surrounding media causing transients to pass rapidly. %Other considerations not accounted for are the effect of internal texture variables within the deforming granular system and nonlocal effects due to finite grain size \edit{\cite{goldsmith2013drag}}.  

A more complex scaling relationship may be needed to go beyond monotonic driving, e.g. oscillatory wheel motion, where the reaction force of the material on the intruder could be more influenced by transient phenomena.  Within this caveat, the scaling could be extended to other types of \edit{locomotion, such as undulating self-propulsion}, with more moving parts by adding additional non-dimensional groups for each new degree of freedom. Testing these extended scaling laws would be important future work.

%In addition to further exploring this rate effect, we also hope to explore potential optimization routines of wheel shapes to allow for more efficient locomotion.

%The scaling constant $\zeta$ could be further expanded, though, using Plasticity theory.  Within Plasticity, it is understood that the only external properties that affect the force of a body intruding into a granular material is the density of the sand $\rho$, the gravitational acceleration $g$, the coefficient of friction of the sand $\mu$, and the coefficient of sand material interface $\mu_w$.  Therefore, the RFT scaling constant $\zeta$ could be rewritten in these terms.  Using dimensional analysis, the relation shown in equation \ref{zetascal} was derived. 
%
%\begin{equation}
%\zeta=\rho g\hat{\zeta}(\mu,\mu_w)
%\label{zetascal}
%\end{equation}
%
%Therefore, the global scaling relations can be rewritten as shown in equations \ref{scalv2} and \ref{scalp2}.
%
%\begin{equation}
%v_{avg}=L\omega\tilde{\psi_v}(f,\frac{g}{L \omega^2},\frac{mg}{F_{draw}},\frac{\rho\hat{\zeta}(\mu,\mu_w) D L^2}{m})
%\label{scalv2}
%\end{equation}
%
%\begin{equation}
%P_{avg}=mL^2\omega^3\tilde{\psi_p}(f,\frac{g}{L \omega^2},\frac{mg}{F_{draw}},\frac{\rho\hat{\zeta}(\mu,\mu_w) D L^2}{m})
%\label{scalp2}
%\end{equation}
%
%\subsection{Functional Scaling Law}
%\indent \indent Using the scaling relations defined in equations \ref{scalv2} and \ref{scalp2} a functional scaling law can be found.  First, to reduce the complexity of the law, the shape and type of sand are assumed to be constant.  This means the wheel shape parameter ($f$) and the sand type ($\hat{\zeta}(\mu,\mu_w)$ and $\rho$) remain constant, resulting in a three degree of freedom system.  We therefore introduce three scaling factors $A$, $B$, and $C$ as defined in equation \ref{scale factors}.
%
%\begin{equation}
%\bar{g}=Ag, \bar{L}=BL, \bar{m}=Cm
%\label{scale factors}
%\end{equation}
%
%Utilizing the global scaling relations, we can solve for the scale factors for the other non-constant inputs as shown in equation \ref{scale fac2}.
%
%\begin{equation}
%\bar{\omega}=\sqrt{\frac{A}{B}}\omega, \bar{F_{draw}}=CAF_{draw}, \bar{D}=\frac{C}{B^2}D
%\label{scale fac2}
%\end{equation}
%
%Inputting these scale factors back into the global scaling relations, we can derive the functional scaling laws as shown in equations \ref{scalelawv} and \ref{scalelawp}.
%
%\begin{equation}
%\begin{multlined}
%v_{avg}(\omega,L,f,m,g,F_{draw},\hat{\zeta}(\mu,\mu_w),D,\rho)= \\
%\frac{1}{\sqrt{AB}}v_{avg}(\sqrt{\frac{A}{B}}\omega,BL,f,Cm,Ag,CAF_{draw},\hat{\zeta}(\mu,\mu_w),\frac{C}{B^2}D,\rho)
%\label{scalelawv}
%\end{multlined}
%\end{equation}
%
%\begin{equation}
%\begin{multlined}
%P_{avg}(\omega,L,f,m,g,F_{draw},\hat{\zeta}(\mu,\mu_w),D,\rho)= \\
%\frac{1}{CA\sqrt{AB}}P_{avg}(\sqrt{\frac{A}{B}}\omega,BL,f,Cm,Ag,CAF_{draw},\hat{\zeta}(\mu,\mu_w),\frac{C}{B^2}D,\rho)
%\label{scalelawp}
%\end{multlined}
%\end{equation}
%
%This means for instance, to test how a wheel design would behave on Mars, one could test a different design on Earth as long as the scaling law is fulfilled.  So, to simulate changes in the gravitational acceleration, one needs only to change accordingly the other inputs $\omega$, $L$, $m$, $F_{draw}$, and $D$.
%
%\subsection{Non-Intuitiveness of Scaling Law}
%\indent \indent As a first check of the scaling laws derived in equations \ref{scalelawv} and \ref{scalelawp}, the case where a cylindrical wheel is rolling without slipping on a solid surface ($\zeta\rightarrow\infty$) can be used.  In this case, one would expect the velocity of both cases to scale purely with the length over the time scale.  Therefore, multiplying the length $BL$ by the rotational velocity $\sqrt{\frac{A}{B}}\omega$ gives an expected scaling of $\sqrt{AB}$, which is what is found.  The power, in the rolling without slipping case, would be $0$ for both wheels, so the scaling law again holds.
%\newline
%\indent Similarly, in the case of a non-round wheel rolling without slipping on a solid surface, one would expect that the velocity would scale the same as the cylindrical wheel.  The power, though, would be expected to scale with the torque, defined as the length $BL$ multiplied by the mass $Cm$ and gravity $Ag$, multiplied by the rotational velocity $\sqrt{\frac{A}{B}}\omega$.  Therefore, the expected power scaling would be $CA\sqrt{AB}$, which is exactly what is found.
%\newline
%\indent Although these scaling relations are valid in both these cases, they are not intuitive.  For instance, one would expect the velocity to scale with $\omega$ and $L$, but it is not expected that changing the mass, gravity, and thickness of the wheel would have any effect.  This relationship between the mass $m$ and the thickness $D$ is especially interesting, because it seems to come directly from RFT.  This scaling could not be found using only Plasticity.
%
%This scaling effect can not only be applied to the average output velocity and power, but it can also be used to find a parameter at any instant of time.  This means, with the velocity, position, and time at every step of two scalable designs, one can alter the time and multiply by a scaling factor to replicate the entire time dependent trajectory.  From the dimensional analysis, we know that the time scales with the rotational velocity $\omega$, such that the same scale factor applied to $\omega$ can be applied to the time $t$.  This effect is the same for the position, which is dimensionally similar to the length $L$ and scales as such.  The velocity as well scales with the length component divided by the time.  Therefore, the position and velocity trajectories can be written as shown in equations \ref{timescalepos} and \ref{timescalevel}, where the two wheels are defined by the subscripts $m$ and $e$, which refers to the Mars and Earth scaling example in the previous section. 
%
%\begin{equation}
%x_m(t_m)=\frac{1}{B}x_e(\sqrt{\frac{A}{B}}t_e)
%\label{timescalepos}
%\end{equation}
%
%\begin{equation}
%v_m(t_m)=\frac{1}{\sqrt{AB}}v_e(\sqrt{\frac{A}{B}}t_e)
%\label{timescalevel}
%\end{equation}
%
%
%
%
%
%\begin{figure}
%\includegraphics[width=0.3\textwidth]{Pictures/4BarWheels.jpg}
%\caption{Four Bar Wheels}
%\label{figFBW}
%\end{figure}
%
%
%\begin{figure}
%\includegraphics[width=0.49\textwidth]{Pictures/AvgSlapMassPower.jpg}
%\caption{Average Slap Mass Power}
%\label{fig05}
%\end{figure}
%
%\begin{figure}
%\includegraphics[width=0.49\textwidth]{Pictures/AvgSlapMassVel.jpg}
%\caption{Average Slap Mass Velocity}
%\label{fig06}
%\end{figure}
%
%\begin{figure}
%\includegraphics[width=0.49\textwidth]{Pictures/AvgSlapGravPower.jpg}
%\caption{Average Slap Grav Power}
%\label{fig07}
%\end{figure}
%
%\begin{figure}
%\includegraphics[width=0.49\textwidth]{Pictures/AvgSlapGravVel.jpg}
%\caption{Average Slap Grav Velocity}
%\label{fig08}
%\end{figure}





%\begin{thebibliography}{20}
%\newpage
\bibliography{granular}
%\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
%\expandafter\ifx\csname bibnamefont\endcsname\relax
%  \def\bibnamefont#1{#1}\fi
%\expandafter\ifx\csname bibfnamefont\endcsname\relax
%  \def\bibfnamefont#1{#1}\fi
%\expandafter\ifx\csname citenamefont\endcsname\relax
%  \def\citenamefont#1{#1}\fi
%\expandafter\ifx\csname url\endcsname\relax
%  \def\url#1{\texttt{#1}}\fi
%\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
%\providecommand{\bibinfo}[2]{#2}
%\providecommand{\eprint}[2][]{\url{#2}}

%\bibitem{ahlvin}
%Ahlvin, Richard B., and Peter W. Haley. Nato Reference Mobility Model: Edition II. NRMM User's Guide. US Army Engineer Waterways Experiment Station, 1992.
%
%\bibitem{Ding}
%Ding, Yang, Nick Gravish, and Daniel I. Goldman. "Drag induced lift in granular media." Physical Review Letters 106.2 (2011): 028001.
%
%\bibitem{Goldmanpropul}
%Goldman, Daniel I. "Colloquium: Biophysical principles of undulatory self-propulsion in granular media." Reviews of Modern Physics 86.3 (2014): 943.
%
%\bibitem{Hennan}
%Henann, David L., and Ken Kamrin. "A predictive, size-dependent continuum model for dense granular flows." Proceedings of the National Academy of Sciences 110.17 (2013): 6730-6735.
%
%\bibitem{Hill}
%Hill, Rodney. The mathematical theory of plasticity. Vol. 11. Oxford university press, 1998.
%
%\bibitem{Kamrin}
%Kamrin, Ken. "Nonlinear elasto-plastic model for dense granular flow." International Journal of Plasticity 26.2 (2010): 167-188.
%
%\bibitem{Koval}
%Kamrin, Ken, and Georg Koval. "Nonlocal constitutive relation for steady granular flow." Physical Review Letters 108.17 (2012): 178301.
%
%\bibitem{}
%
%\bibitem{}
%
%\bibitem{}
%
%\bibitem{}
%
%\bibitem{supp}
%See Supplemental Material
%%See Supplemental Material for a detailed discussion of the intruder boundary conditions, a justification of our choice of $A$, a mesh sensitivity study, additional simulated creep data, and explanations for exponential creep-rate behavior.

%\end{thebibliography}
%
%\newpage
%\textbf{Appendix}
%\kcom{Have a look at how I squeezed these into the main body and modify as you see fit.  After that we may not need these appendices.}
%\jcom{I think it works a lot better the way you squeezed it into the paper, so I don't see any reason we need the appendix.  I made the comments above about time being left out though in case you did that accidentally.}
%
%\begin{multline}
%P_s\left(t_s,f_s,M_s,L_s,\omega_s,D_s\right)=\\
%\frac{1}{C\sqrt{B}}P_b\left(\sqrt{B}t_s,f_s,CM_s,BL_s,\frac{\omega_s}{\sqrt{B}},\frac{CD_s}{B^2}\right)
%\label{FuncScalP1}
%\end{multline}
%\begin{multline}
%V_s\left(t_s,f_s,M_s,L_s,\omega_s,D_s\right)=\\
%\frac{1}{\sqrt{B}}V_b\left(\sqrt{B}t_s,f_s,CM_s,BL_s,\frac{\omega_s}{\sqrt{B}},\frac{CD_s}{B^2}\right)
%\label{FuncScalV1}
%\end{multline}
%\jcom{This one is the functional scaling form without gravity scaling as one equation.  AKA the first two equations combined into one.}
%\begin{multline}
%\left[P_s,V_s\right]\left(t_s,f_s,M_s,L_s,\omega_s,D_s\right)=\\
%\left[\frac{P_b}{C\sqrt{B}},\frac{V_b}{\sqrt{B}}\right]\left(\sqrt{B}t_s,f_s,CM_s,BL_s,\frac{\omega_s}{\sqrt{B}},\frac{CD_s}{B^2}\right)
%\label{FuncScalP1}
%\end{multline}
%
%\jcom{These two are the functional scaling form WITH gravity scaling as two equations.}
%\begin{multline}
%V_s\left(t_s,f_s,M_s,L_s,\omega_s,g_s,\hat\xi_s,D_s,\rho_s\right)=\frac{V_b}{\sqrt{AB}}\\
%\left(\sqrt{\frac{B}{A}}t_s,f_s,CM_s,BL_s,\sqrt{\frac{A}{B}}\omega_s,Ag_s,\hat\xi_s,\frac{C}{B^2}D_s,\rho_s\right)
%\label{FuncScalV1}
%\end{multline}
%\begin{multline}
%P_s\left(t_s,f_s,M_s,L_s,\omega_s,g_s,\hat\xi_s,D_s,\rho_s\right)=\frac{P_b}{CA\sqrt{AB}}\\
%\left(\sqrt{\frac{B}{A}}t_s,f_s,CM_s,BL_s,\sqrt{\frac{A}{B}}\omega_s,Ag_s,\hat\xi_s,\frac{C}{B^2}D_s,\rho_s\right)
%\label{FuncScalV1}
%\end{multline}
%
%\jcom{This one is the functional scaling form WITH gravity scaling as one equation.  AKA the previous two equations combined into one.}
%\begin{multline}
%\left[P_s,V_s\right]\left(t_s,f_s,M_s,L_s,\omega_s,g_s,\hat\xi_s,D_s,\rho_s\right)=\left[\frac{P_b}{CA\sqrt{AB}},\frac{V_b}{\sqrt{AB}}\right]\\\left(\sqrt{\frac{B}{A}}t_s,f_s,CM_s,BL_s,\sqrt{\frac{A}{B}}\omega_s,Ag_s,\hat\xi_s,\frac{C}{B^2}D_s,\rho_s\right)
%\label{FuncScalP1}
%\end{multline}
%

\end{document}

RFT_ReSubmit3.bbl
\begin{thebibliography}{30}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\expandafter\ifx\csname bibnamefont\endcsname\relax
  \def\bibnamefont#1{#1}\fi
\expandafter\ifx\csname bibfnamefont\endcsname\relax
  \def\bibfnamefont#1{#1}\fi
\expandafter\ifx\csname citenamefont\endcsname\relax
  \def\citenamefont#1{#1}\fi
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
\providecommand{\bibinfo}[2]{#2}
\providecommand{\eprint}[2][]{\url{#2}}

\bibitem[{\citenamefont{Goddard}(2014)}]{goddard2014continuum}
\bibinfo{author}{\bibfnamefont{J.~D.} \bibnamefont{Goddard}},
  \bibinfo{journal}{Applied Mechanics Reviews} \textbf{\bibinfo{volume}{66}},
  \bibinfo{pages}{050801} (\bibinfo{year}{2014}).

\bibitem[{\citenamefont{Schofield and Wroth}(1968)}]{schofield1968critical}
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Schofield}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Wroth}}
  (\bibinfo{year}{1968}).

\bibitem[{\citenamefont{Wood}(1990)}]{wood1990soil}
\bibinfo{author}{\bibfnamefont{D.~M.} \bibnamefont{Wood}},
  \emph{\bibinfo{title}{Soil behaviour and critical state soil mechanics}}
  (\bibinfo{publisher}{Cambridge university press}, \bibinfo{year}{1990}).

\bibitem[{\citenamefont{Roux and Radjai}(1998)}]{roux1998texture}
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Roux}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Radjai}}, in
  \emph{\bibinfo{booktitle}{Physics of dry granular media}}
  (\bibinfo{publisher}{Springer}, \bibinfo{year}{1998}), pp.
  \bibinfo{pages}{229--236}.

\bibitem[{\citenamefont{Depken et~al.}(2007)\citenamefont{Depken, Lechman, van
  Hecke, van Saarloos, and Grest}}]{depken07}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Depken}},
  \bibinfo{author}{\bibfnamefont{J.~B.} \bibnamefont{Lechman}},
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{van Hecke}},
  \bibinfo{author}{\bibfnamefont{W.}~\bibnamefont{van Saarloos}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{G.~S.} \bibnamefont{Grest}},
  \bibinfo{journal}{Europhys. Lett.} \textbf{\bibinfo{volume}{78}},
  \bibinfo{pages}{58001} (\bibinfo{year}{2007}).

\bibitem[{\citenamefont{Weinhart et~al.}(2013)\citenamefont{Weinhart, Hartkamp,
  Thornton, and Luding}}]{weinhart2013coarse}
\bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Weinhart}},
  \bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Hartkamp}},
  \bibinfo{author}{\bibfnamefont{A.~R.} \bibnamefont{Thornton}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Luding}},
  \bibinfo{journal}{Physics of Fluids (1994-present)}
  \textbf{\bibinfo{volume}{25}}, \bibinfo{pages}{070605}
  (\bibinfo{year}{2013}).

\bibitem[{\citenamefont{S}(2008)}]{luding08}
\bibinfo{author}{\bibfnamefont{L.}~\bibnamefont{S}},
  \bibinfo{journal}{Particul. Sci. Technol.} \textbf{\bibinfo{volume}{26}},
  \bibinfo{pages}{33} (\bibinfo{year}{2008}).

\bibitem[{\citenamefont{Bagnold}(1954)}]{bagnold54}
\bibinfo{author}{\bibfnamefont{R.~A.} \bibnamefont{Bagnold}},
  \bibinfo{journal}{Proc. Roy. Soc. London Ser. A}
  \textbf{\bibinfo{volume}{225}} (\bibinfo{year}{1954}).

\bibitem[{\citenamefont{da~Cruz et~al.}(2005)\citenamefont{da~Cruz, Emam,
  Prochnow, Roux, and Chevoir}}]{dacruz05}
\bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{da~Cruz}},
  \bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Emam}},
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Prochnow}},
  \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Roux}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Chevoir}},
  \bibinfo{journal}{Phys. Rev. E.} \textbf{\bibinfo{volume}{72}},
  \bibinfo{pages}{021309} (\bibinfo{year}{2005}).

\bibitem[{\citenamefont{MiDi}(2004)}]{midi2004dense}
\bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{MiDi}}, \bibinfo{journal}{The
  European Physical Journal E} \textbf{\bibinfo{volume}{14}},
  \bibinfo{pages}{341} (\bibinfo{year}{2004}).

\bibitem[{\citenamefont{Jop et~al.}(2006)\citenamefont{Jop, Forterre, and
  Pouliquen}}]{jop06}
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Jop}},
  \bibinfo{author}{\bibfnamefont{Y.}~\bibnamefont{Forterre}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{O.}~\bibnamefont{Pouliquen}},
  \bibinfo{journal}{Nature} \textbf{\bibinfo{volume}{441}},
  \bibinfo{pages}{727} (\bibinfo{year}{2006}).

\bibitem[{\citenamefont{Kamrin and Koval}(2012)}]{kamrin2012nonlocal}
\bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Kamrin}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Koval}},
  \bibinfo{journal}{Physical Review Letters} \textbf{\bibinfo{volume}{108}},
  \bibinfo{pages}{178301} (\bibinfo{year}{2012}).

\bibitem[{\citenamefont{Henann and Kamrin}(2013)}]{henann2013predictive}
\bibinfo{author}{\bibfnamefont{D.~L.} \bibnamefont{Henann}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Kamrin}},
  \bibinfo{journal}{Proceedings of the National Academy of Sciences}
  \textbf{\bibinfo{volume}{110}}, \bibinfo{pages}{6730} (\bibinfo{year}{2013}).

\bibitem[{\citenamefont{Kamrin and Henann}(2015)}]{kamrin2015nonlocal}
\bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Kamrin}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{D.~L.} \bibnamefont{Henann}},
  \bibinfo{journal}{Soft matter} \textbf{\bibinfo{volume}{11}},
  \bibinfo{pages}{179} (\bibinfo{year}{2015}).

\bibitem[{\citenamefont{S.~Mohan et~al.}(2002)\citenamefont{S.~Mohan, K.~Rao,
  and Nott}}]{mohan02}
\bibinfo{author}{\bibfnamefont{L.}~\bibnamefont{S.~Mohan}},
  \bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{K.~Rao}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{P.~R.} \bibnamefont{Nott}},
  \bibinfo{journal}{J. Fluid Mech.} \textbf{\bibinfo{volume}{457}},
  \bibinfo{pages}{377} (\bibinfo{year}{2002}).

\bibitem[{\citenamefont{Aranson and Tsimring}(2002)}]{aranson02}
\bibinfo{author}{\bibfnamefont{I.~S.} \bibnamefont{Aranson}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{L.~S.} \bibnamefont{Tsimring}},
  \bibinfo{journal}{Phys. Rev. E} \textbf{\bibinfo{volume}{65}},
  \bibinfo{pages}{061303} (\bibinfo{year}{2002}).

\bibitem[{\citenamefont{Li et~al.}(2013)\citenamefont{Li, Zhang, and
  Goldman}}]{li13}
\bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Li}},
  \bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Zhang}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{D.~I.} \bibnamefont{Goldman}},
  \bibinfo{journal}{Science} \textbf{\bibinfo{volume}{339}},
  \bibinfo{pages}{1408} (\bibinfo{year}{2013}).

\bibitem[{\citenamefont{Lighthill}(1975)}]{lighthill75}
\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Lighthill}},
  \emph{\bibinfo{title}{Mathematical Biofluiddynamics}}
  (\bibinfo{publisher}{SIAM}, \bibinfo{year}{1975}).

\bibitem[{\citenamefont{Goldman}(2014)}]{goldman14}
\bibinfo{author}{\bibfnamefont{D.~I.} \bibnamefont{Goldman}},
  \bibinfo{journal}{Reviews of Modern Physics} \textbf{\bibinfo{volume}{86}},
  \bibinfo{pages}{943} (\bibinfo{year}{2014}).

\bibitem[{\citenamefont{Ding et~al.}(2011)\citenamefont{Ding, Gravish, and
  Goldman}}]{ding11}
\bibinfo{author}{\bibfnamefont{Y.}~\bibnamefont{Ding}},
  \bibinfo{author}{\bibfnamefont{N.}~\bibnamefont{Gravish}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{D.~I.} \bibnamefont{Goldman}},
  \bibinfo{journal}{Physical review letters} \textbf{\bibinfo{volume}{106}},
  \bibinfo{pages}{028001} (\bibinfo{year}{2011}).

\bibitem[{\citenamefont{Askari and Kamrin}(2016)}]{askari2016intrusion}
\bibinfo{author}{\bibfnamefont{H.}~\bibnamefont{Askari}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Kamrin}},
  \bibinfo{journal}{Nature Materials} \textbf{\bibinfo{volume}{15}},
  \bibinfo{pages}{1274} (\bibinfo{year}{2016}).

\bibitem[{\citenamefont{Wong and Reece}(1967)}]{wong67}
\bibinfo{author}{\bibfnamefont{J.-Y.} \bibnamefont{Wong}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Reece}},
  \bibinfo{journal}{Journal of Terramechanics} \textbf{\bibinfo{volume}{4}},
  \bibinfo{pages}{81} (\bibinfo{year}{1967}).

\bibitem[{\citenamefont{Ahlvin and Haley}(1992)}]{ahlvin92}
\bibinfo{author}{\bibfnamefont{R.~B.} \bibnamefont{Ahlvin}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{P.~W.} \bibnamefont{Haley}},
  \emph{\bibinfo{title}{Nato Reference Mobility Model: Edition II. NRMM User's
  Guide}} (\bibinfo{publisher}{US Army Engineer Waterways Experiment Station},
  \bibinfo{year}{1992}).

\bibitem[{\citenamefont{Bekker}(1956)}]{bekker}
\bibinfo{author}{\bibfnamefont{M.~G.} \bibnamefont{Bekker}},
  \emph{\bibinfo{title}{Theory of Land Locomotion}}
  (\bibinfo{publisher}{University of Michigan Press}, \bibinfo{year}{1956}).

\bibitem[{\citenamefont{Nedderman}(1992)}]{nedderman}
\bibinfo{author}{\bibfnamefont{R.~M.} \bibnamefont{Nedderman}},
  \emph{\bibinfo{title}{Statics and Kinematics of Granular Materials}}
  (\bibinfo{publisher}{Cambridge University Press}, \bibinfo{year}{1992}).

\bibitem[{\citenamefont{Senatore and Iagnemma}(2014)}]{senatore2014analysis}
\bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Senatore}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Iagnemma}},
  \bibinfo{journal}{Journal of Terramechanics} \textbf{\bibinfo{volume}{51}},
  \bibinfo{pages}{1} (\bibinfo{year}{2014}).

\bibitem[{\citenamefont{Plimpton}(1995)}]{plimpton1995fast}
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Plimpton}},
  \bibinfo{journal}{Journal of computational physics}
  \textbf{\bibinfo{volume}{117}}, \bibinfo{pages}{1} (\bibinfo{year}{1995}).

\bibitem[{\citenamefont{Koval et~al.}(2009)\citenamefont{Koval, Roux, Corfdir,
  and Chevoir}}]{koval2009annular}
\bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Koval}},
  \bibinfo{author}{\bibfnamefont{J.-N.} \bibnamefont{Roux}},
  \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Corfdir}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Chevoir}},
  \bibinfo{journal}{Physical Review E} \textbf{\bibinfo{volume}{79}},
  \bibinfo{pages}{021306} (\bibinfo{year}{2009}).

\bibitem[{\citenamefont{Lagr{\'e}e et~al.}(2011)\citenamefont{Lagr{\'e}e,
  Staron, and Popinet}}]{lagree2011granular}
\bibinfo{author}{\bibfnamefont{P.-Y.} \bibnamefont{Lagr{\'e}e}},
  \bibinfo{author}{\bibfnamefont{L.}~\bibnamefont{Staron}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Popinet}},
  \bibinfo{journal}{Journal of Fluid Mechanics} \textbf{\bibinfo{volume}{686}},
  \bibinfo{pages}{378} (\bibinfo{year}{2011}).

\bibitem[{\citenamefont{Goldsmith et~al.}(2013)\citenamefont{Goldsmith, Guo,
  Hunt, Tao, and Koehler}}]{goldsmith2013drag}
\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Goldsmith}},
  \bibinfo{author}{\bibfnamefont{H.}~\bibnamefont{Guo}},
  \bibinfo{author}{\bibfnamefont{S.~N.} \bibnamefont{Hunt}},
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Tao}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Koehler}},
  \bibinfo{journal}{Physical Review E} \textbf{\bibinfo{volume}{88}},
  \bibinfo{pages}{030201} (\bibinfo{year}{2013}).

\end{thebibliography}
