# Covering convection with thermal blankets: formation of supercontinents

Jinzi Mac Huang1,2<EMAIL_ADDRESS>1NYU-ECNU Institute of Physics and
Institute of Mathematical Sciences, New York University Shanghai, Shanghai,
China 2Applied Math Lab, Courant Institute, New York University, New York, USA

###### Abstract

The continental plates of Earth are known to drift over a geophysical
timescale, and their interactions have lead to some of the most spectacular
geoformations of our planet while also causing natural disasters such as
earthquakes and volcanic activity. Understanding the dynamics of interacting
continental plates is thus significant. In this work, we present a fluid
mechanical investigation of the plate motion, interaction, and dynamics.
Through numerical experiments, we examine the coupling between a convective
fluid and plates floating on top of it. With physical modeling, we show the
coupling is both mechanical and thermal, leading to the thermal blanket
effect: the floating plate is not only transported by the fluid flow beneath,
it also prevents the heat from leaving the fluid, leading to a convective flow
that further affects the plate motion. By adding several plates to such a
coupled fluid-structure interaction, we also investigate how floating plates
interact with each other and show that, under proper conditions, small plates
can converge into a supercontinent.

###### keywords:

Authors should not enter keywords on the manuscript, as these must be chosen
by the author during the online submission process and will then be added
during the typesetting process (see
http://journals.cambridge.org/data/relatedlink/jfm-keywords.pdf for the full
list)

##  1 Introduction

Figure 1:  Rayleigh-Bénard convection coupled to a free-floating plate leads
to different dynamics of plate motion. (a) Schematics of the interaction
between Rayleigh-Bénard convection and the floating plate. The fluid is heated
from below and has an open free surface, the floating plate on this free
surface is transported by the fluid force exerted from below. (b) Different
states of plate motion. A small plate with $d/D=0.2$ oscillates between two
sidewalls of the convection cell, while a big plate with $d/D=0.7$ is trapped
in the middle of the convection cell. Here $L=(D-d)/2$ is the bound of plate
center $x_{p}$. (c) Flow visualization reveals two counter-rotating large-
scale circulations when the plate is located at the center of the convection
cell. Image credit: Zhong & Zhang (2007b); Huang et al. (2018).

Fluid-structure interactions appear at many different scales on our planet,
and perhaps the largest one is the continental plate tectonics (Plummer et
al., 2001). It is believed that this tectonic motion results from the thermal
convection in Earth’s mantle (Kious & Tilling, 1996), where the mantle
materials behave like a fluid (Turcotte & Schubert, 2002) that is heated from
the core and cooled at the surface of Earth. Under gravity, this configuration
of heating and cooling leads to thermal convection, whose circulation provides
forcing to the continental plates through shearing. This is considered to be
the simplest picture of plate tectonics, however many details of the plate
dynamics, like the periodic formation of supercontinents and the associated
geological Wilson cycle (Burke, 2011), require further investigation.

Laboratory experiments have proven to be an effective way of understanding the
fluid-structure interaction behind the plate tectonics, with many successful
studies that couple the thermally convective fluid to solid structures (Elder,
1967; Howard et al., 1970; Whitehead, 1972). Aimed at recovering the plate
dynamics and understanding the associated fluid-structure interactions, a
series of laboratory experiments was later conducted by Zhang & Libchaber
(2000); Zhong & Zhang (2005, 2007a, 2007b). Shown in Fig. 1(a), this
experiment employs water as the working fluid, and a heater beneath provides
heating while the ventilation at water-air interface provides cooling,
resulting in Rayleigh-Bénard convection (Ahlers et al., 2009). A floating
plate of size $d$ is carefully placed on top of this convective domain of
total length $D$. This moving plate with center location $x_{p}$ has only
fluid force acting on it, unless it hits the wall on the left or the right.
Large-scale circulations (Araujo et al., 2005; Brown & Ahlers, 2007; Moore &
Huang, 2023) are observed to form in the convective fluid, as shown in Fig.
1(a) and Fig. 1(c). Depending on the location $x_{p}$, the plate can be either
transported by the circulating fluid, or situated on top of a converging or a
diverging center of the surface flow [this is the case shown in Fig. 1(a)].
Interestingly, different plate motions are observed depending on its size:
When the ratio between plate size $d$ and tank width $D$ is smaller than a
critical number near $0.58$, the plate oscillates between the two sidewalls as
shown in Fig. 1(b); When this ratio is above the critical value, the plate is
trapped at the center of tank, as shown in Fig. 1(b) and Fig. 1(c).

Zhong, Zhang, and others investigated this transition of behaviors, and they
discovered that the so-called thermal blanket effect is responsible here
(Zhang & Libchaber, 2000; Zhong & Zhang, 2005, 2007a, 2007b; Huang et al.,
2018; Mao et al., 2019; Mao, 2021; Lowenstein, 2024). In this theory, the
floating plate serves as an insulator (like a blanket) on top of the
convective fluid, hence the fluid beneath becomes warmer due to the lack of
ventilation. The warm fluid then rises, creating a diverging surface flow as
shown in Fig. 1(a) that can transport the plate. The coupling between the
fluid and the floating plate therefore goes two ways: the plate modifies the
flow temperature and leads to thermal convection; the formed convective flows
transport the floating plate. Their interplay leads to nontrivial dynamics of
the plate shown in Fig. 1(b), and the physically inspired Zhong-Zhang model
(Zhong & Zhang, 2005, 2007a, 2007b) successfully captures the transition of
dynamics. Recently, more careful investigations on the Zhong-Zhang model have
lead to new advancements in the stochastic (Huang et al., 2018) and dynamical
(Lowenstein, 2024) modeling of fluid-structure interactions.

While the laboratory experiments are conducted in a domain of fluid with
finite size, numerical simulations can be conducted in a domain that resembles
the mantle of Earth. The numerical work of Gurnis (1988) provides one of the
first time-dependent simulations of continental drift, where the fluid domain
is 2-dimensional and periodic. After this, many more numerical works have
investigated the details of continental drift (Zhong & Gurnis, 1993; Lowman &
Jarvis, 1993, 1995; Lowman & Gable, 1999; Lowman & Jarvis, 1999; Lowman et
al., 2001; Zhong et al., 2000), engaging higher resolutions, more detailed
modeling of fluid-structure interactions, and 3-dimensional simulations of the
interior of Earth. In recent years, the mobility of the continental plate has
become a focus of numerical study, where persistent motion is observed for
larger plates while small plates tend to move sporadically (Gurnis, 1988;
Whitehead & Behn, 2015; Mao et al., 2019; Mao, 2021). In these works, the
thermal blanket effect is once again recognized as an important factor causing
the diverse plate dynamics.

This work is a continuation of an earlier investigation, Huang (2024), where
the thermal and mechanical coupling between one floating plate and convective
fluid is modelled through a simple stochastic model. This model shows that the
covering ratio Cr, defined as the ratio of the plate area to the total surface
area, is a direct measure of the thermal blanket effect. A critical covering
ratio $\mbox{{Cr}}^{*}$ is identified to distinguish the dynamics of the
plate: For a small plate with $\mbox{{Cr}}<\mbox{{Cr}}^{*}$, the plate is
passive to the flow field and exhibits little motion; For a plate with
$\mbox{{Cr}}>\mbox{{Cr}}^{*}$, the strong thermal blanket effect leads to
persistent plate translation.

In this work, we first introduce an efficient 2-dimensional spectral solver
that can accurately capture the motions and interactions of multiple floating
plates on top of Rayleigh-Bénard convection. In a periodic domain shown in
Fig. 2, this solver can handle the Navier-Stokes flows presented in laboratory
experiments, with simple modifications available for the geophysical Stokes
flows in the mantle. Moreover, multiple floating plates can be simulated as
fast as a single plate problem, as the floating plates are simply treated as
an area with different boundary conditions. A specially-tailored spectral
solver handles the resulting inhomogeneous Robin conditions for both the
temperature and the stream function, allowing for efficient time-stepping and
spectral accuracy. This enables us to systematically introduce 1, 2, and many
floating plates, and to show how the thermal blanket effect dictates their
interactions with the convective flow beneath and each other. The covering
ratio Cr is once again identified as a key factor affecting the plate dynamics
and the stable formation of supercontinents.

This paper is arranged as follows: In Section 2, we will mathematically
formulate the Rayleigh-Bénard convection and its coupling to the plate motion;
In Section 3, a numerical scheme and its implementation for solving this free-
boundary problem will be introduced; In Section 4, numerical simulations of
single, double, and multiple plate dynamics will be included and discussed;
Finally, we will discuss extensions and future works in Section 5.

##  2 Mathematical formulation

###  2.1 Flow and temperature equations

We consider a dimensionless set of equations by rescaling the length scale by
the cell height $H$, the time scale by the diffusion time $H^{2}/\kappa$
($\kappa$ is thermal diffusivity), and temperature by the temperature
difference $\Updelta T$ between the heater and free surface. The $x$ direction
of the fluid domain is periodic with period $\Gamma=D/H$ ($D$ is the domain
width), so the overall computational domain is $x\in(0,\Gamma)$, $y\in(0,1)$
as shown in Fig. 2. With the Boussinesq approximation, the resulting PDEs for
flow speed $\bm{u}=(u,v)$, pressure $p$, and temperature $\theta\in[0,1]$ are

$\displaystyle\frac{D\bm{u}}{Dt}=-\nabla
p+\text{Pr}{}\,\nabla^{2}\bm{u}+\text{Ra}{}\,\text{Pr}{}\,\theta,$ (1)
$\displaystyle\nabla\cdot\bm{u}=0,$ (2)
$\displaystyle\frac{D\theta}{Dt}=\nabla^{2}\theta.$ (3)

Two dimensionless numbers appear during this non-dimensionalization: the
Rayleigh number $\text{Ra}{}=\alpha g\Updelta TH^{3}/\nu\kappa$ and the
Prandtl number $\text{Pr}{}=\nu/\kappa$, with $\nu,\alpha$ and $g$ denoting
kinematic viscosity, thermal expansion coefficient of the fluid, and the
acceleration due to gravity. We have further assumed that the physical
properties of the fluid depend on temperature weakly, so Ra and Pr do not
depend on $\theta$.

As our simulation is 2-dimensional, it is convenient to write the Navier-
Stokes equation in a vorticity & steam function format (Peyret, 2002),

$\displaystyle\frac{D\omega}{Dt}$
$\displaystyle=\text{Pr}{}\nabla^{2}\omega+\text{Pr}{}\text{Ra}{}\frac{\partial\theta}{\partial
x},$ (4) $\displaystyle-\nabla^{2}\psi$ $\displaystyle=\omega,\ \
\bm{u}=\nabla_{\perp}\psi,$ (5) $\displaystyle\frac{D\theta}{Dt}$
$\displaystyle=\nabla^{2}\theta.$ (6)

Instead of directly solving for $\bm{u}$ and $p$, the $z$-component of
vorticity $\omega=\mathbf{k}\cdot\nabla\times\bm{u}$ and the stream function
defined by $\bm{u}=\nabla_{\perp}\psi=(\psi_{y},-\psi_{x})$ are solved first.

###  2.2 Boundary conditions

Figure 2: Schematics of the floating plate problem. The fluid domain $\Omega$
is heated from the bottom surface $\partial\Omega_{0}$ and has an open surface
on top ($\partial\Omega_{1}$), floating plates $P_{1},P_{1},P_{2},\dots$ cover
part of this open surface and shield the heat from leaving the fluid.

While Eqs. 4, 5 and 6 are standard for modeling Rayleigh-Bénard convection,
the boundary conditions become complicated when floating plates are present.
Shown in Fig. 2, the fluid domain $\Omega$ is bounded between the bottom
heating wall $y=0$ ($\partial\Omega_{0}$) and the top free surface $y=1$
($\partial\Omega_{1}$). The segments of top surface covered by the floating
plates are labelled as $P_{1},P_{2},\cdots$, whose centers are
$x_{p}^{(1)},x_{p}^{(2)},\cdots$.

The boundary conditions for the bottom heating wall are set straightforwardly
as constant temperature and no-slip,

$\theta=1,\,u=v=0\quad\mbox{at }y=0.$ (7)

For the vorticity-stream function format,

$\theta=1,\,\psi=\psi_{y}=0\quad\mbox{at }y=0.$ (8)

The top condition consists of two types of regions: for the free surface (not
covered by $P_{i}$), the temperature is low and the flow is shear free; for
the region covered by the plate $P_{i}$, the heat flux is zero and the flow
shares the same velocity as the plate. The zero-flux condition originates from
the “thermal blanket” effect caused by the low heat conduction of solids. The
boundary conditions at $y=1$ can be summarized as

$\displaystyle\theta=0,\,u_{y}=v=0\quad$ $\displaystyle\mbox{for }y=1\mbox{
and }x\notin\cup P_{i},$ (9)
$\displaystyle\theta_{y}=0,\,u=u_{p}^{(i)},v=0\quad$ $\displaystyle\mbox{for
}y=1\mbox{ and }x\in P_{i}.$ (10)

Here $u_{p}^{(i)}=\dot{x}_{p}^{(i)}$ is the velocity of $i$th plate $P_{i}$.
For convenience, we can also write the top boundary conditions in a more
compact way,

$\begin{cases}(1-a)\,\theta+a\,\theta_{y}=0\\\
a\,u+(1-a)\,u_{y}=g\quad\quad\mbox{at }y=1.\\\ v=0\end{cases}$ (11)

For the vorticity-stream function format,

$\begin{cases}(1-a)\,\theta+a\,\theta_{y}=0\\\
a\,\psi_{y}+(1-a)\,\psi_{yy}=g\quad\quad\mbox{at }y=1.\\\ \psi=0\end{cases}$
(12)

Here $a(x)=\sum_{i}\mathbb{1}_{x\in P_{i}}$,
$g(x)=\sum_{i}u_{p}^{(i)}\mathbb{1}_{x\in P_{i}}$, and $\mathbb{1}_{x\in
P_{i}}$ is the characteristic function such that

$\mathbb{1}_{x\in P_{i}}=\begin{cases}1\quad\mbox{ if }x\in P_{i},\\\
0\quad\mbox{ otherwise. }\end{cases}$ (13)

On plate $P_{i}$, two types of forces drive its motion: the fluid force
$f^{(i)}$ due to the shear from convective flows; the interacting force
$f_{l}^{(i)}$ or $f_{r}^{(i)}$ when the left or right neighboring plate
($P_{i-1}$ or $P_{i+1}$) makes contact with plate $P_{i}$.

For the fluid force, we simply integrate the fluid shear stress,

$f^{(i)}=-\text{Pr}\int_{P_{i}}u_{y}(x,1,t)\,dx=-\text{Pr}\int_{0}^{\Gamma}u_{y}(x,1,t)\mathbb{1}_{x\in
P_{i}}d\,x.$ (14)

The interaction forces $f_{l}^{(i)}$ and $f_{r}^{(i)}$ are modelled as a
short-range interaction force to ensure a fully-elastic collision between
plates. The numerical implementations will be included in Section 3.6.

Finally, we add all the forces together and evolve the plate location as

$\displaystyle\dot{x}_{p}^{(i)}$ $\displaystyle=u_{p}^{(i)}{},$ (15)
$\displaystyle\dot{u}_{p}^{(i)}$
$\displaystyle=a_{p}^{(i)}=m^{-1}\left[\,f_{l}^{(i)}+f_{r}^{(i)}-\text{Pr}\int_{0}^{\Gamma}u_{y}(x,1,t)\mathbb{1}_{x\in
P_{i}}\,dx\right].$ (16)

Here $a_{p}^{(i)}$ is the acceleration of $P_{i}$ and $m$ is the dimensionless
mass of plate.

##  3 Numerical methods

###  3.1 Time stepping

We discretize time with the second order Adam-Bashforth Backward
Differentiation method (ABBD2). At time step $t_{n}=n\Updelta T$, Eqs. 4, 6
and 5 become

$\displaystyle\nabla^{2}\omega_{n}-\sigma_{1}\omega_{n}$
$\displaystyle=f_{n},$ (17)
$\displaystyle\nabla^{2}\theta_{n}-\sigma_{2}\theta_{n}$
$\displaystyle=h_{n},$ (18) $\displaystyle-\nabla^{2}\psi_{n}$
$\displaystyle=\omega_{n},$ (19)

where

$\displaystyle\sigma_{1}$ $\displaystyle=\frac{3}{2\,\text{Pr}\,\Updelta
t},\quad\sigma_{2}=\frac{3}{2\Updelta t},$ (20) $\displaystyle f_{n}$
$\displaystyle=\text{Pr}^{-1}\left[2(\bm{u}\cdot\nabla\omega)_{n-1}-(\bm{u}\cdot\nabla\omega)_{n-2}\right]$
(21) $\displaystyle\hskip 40.0pt-(2\,\text{Pr}\,\Updelta
t)^{-1}\left(4\omega_{n-1}-\omega_{n-2}\right)-\text{Ra}\,\left(\frac{\partial\theta}{\partial
x}\right)_{n},$ $\displaystyle h_{n}$
$\displaystyle=\left[2(\bm{u}\cdot\nabla\theta)_{n-1}-(\bm{u}\cdot\nabla\theta)_{n-2}\right]-(2\Updelta
t)^{-1}\left(4\theta_{n-1}-\theta_{n-2}\right).$ (22)

Equations 17, 18 and 19 are Helmholtz equations that can be solved by standard
numerical methods (Peyret, 2002), and this implicit-explicit operator
splitting scheme has been used in many numerical studies of thermal convection
(Huang & Zhang, 2022; Huang, 2024). However, modifications have to be made to
accommodate for the inhomogeneous Robin boundary conditions Eq. 12. We will
detail this Helmholtz solver in the next few sections.

In Eqs. 21 and 22, nonlinear terms like $\bm{u}\cdot\nabla\theta$ and
$\bm{u}\cdot\nabla\omega$ are computed pseudo-spectrally, with a simple and
efficient anti-aliasing filter (Hou & Li, 2007). With given initial and
boundary data, (18) can be solved first to obtain $\theta_{n}$, which is
inserted in $f_{n}$ so (17) can be solved next. Finally, (19) is solved with
the known $\omega_{n}$. During time stepping, we typically set $\Updelta
t=5\times 10^{-4}\,\text{Ra}^{-1/2}$ considering that the flow speed
$|\bm{u}|\sim\sqrt{\text{Ra}}$.

After solving for the flow and temperature fields, the acceleration
$a_{p,\,n}^{(i)}$ of plate $P_{i}$ can be computed via Eq. 16, and Eqs. 15 and
16 are integrated with a $2^{\textnormal{nd}}$ order Adam-Bashforth method,

$\displaystyle x_{p,\,n+1}^{(i)}$
$\displaystyle=x_{p,\,n}^{(i)}+\frac{\Updelta
t}{2}\left[3u_{p,\,n}^{(i)}-u_{p,\,n-1}^{(i)}\right],$ (23) $\displaystyle
u_{p,\,n+1}^{(i)}$ $\displaystyle=u_{p,\,n}^{(i)}+\frac{\Updelta
t}{2}\left[3a_{p,\,n}^{(i)}-a_{p,\,n-1}^{(i)}\right].$ (24)

###  3.2 Helmholtz solver for Dirichlet problems

Before constructing the numerical solver for Robin boundary conditions, we
first outline the numerical solver for the following Helmholtz problem,

$\displaystyle\nabla^{2}u-\sigma u$ $\displaystyle=h(x,y),$ (25)
$\displaystyle u(x,0)$ $\displaystyle=g_{0}(x),$ (26) $\displaystyle u(x,1)$
$\displaystyle=g_{1}(x).$ (27)

As the domain is periodic in $x\in(0,\Gamma]$, we can discretize $x$ as $L$
equally-spaced nodes so $x_{l}=l\Updelta x$ where $\Updelta x=\Gamma/L$ and
$l=1,2,\dots L$. We further require $L$ to be odd. We can now approximate the
solution $u(x,y)$ as a truncated Fourier series in $x$,

$u(x_{l},y)=\sum_{k=-\frac{L-1}{2}}^{\frac{L-1}{2}}\hat{u}_{k}(y)\exp{(2\pi\mathbb{i}kl/L)},$
(28)

where $\mathbb{i}$ is the imaginary unit. The Fourier coefficients
$\hat{u}_{k}$ then satisfy the ODE,

$\displaystyle\frac{d^{2}\hat{u}_{k}}{dy^{2}}-\sigma_{k}\hat{u}_{k}$
$\displaystyle=\hat{h}_{k}(y),$ (29) $\displaystyle\hat{u}_{k}(0)$
$\displaystyle=\hat{g}_{0,k},$ (30) $\displaystyle\hat{u}_{k}(1)$
$\displaystyle=\hat{g}_{1,k},$ (31)

where

$\displaystyle\sigma_{k}$
$\displaystyle=\frac{4\pi^{2}k^{2}}{\Gamma^{2}}+\sigma,$ (32)
$\displaystyle\hat{h}_{k}(y)$
$\displaystyle=\frac{1}{L}\sum_{l=1}^{L}h(x_{l},y)\exp{(-2\pi\mathbb{i}kl/L)},$
(33) $\displaystyle\hat{g}_{0,k}$
$\displaystyle=\frac{1}{L}\sum_{l=1}^{L}g_{0}(x_{l})\exp{(-2\pi\mathbb{i}kl/L)},$
(34) $\displaystyle\hat{g}_{1,k}$
$\displaystyle=\frac{1}{L}\sum_{l=1}^{L}g_{1}(x_{l})\exp{(-2\pi\mathbb{i}kl/L)}.$
(35)

The computation of Eqs. 33, 34 and 35 can be done efficiently with the help of
the Fast Fourier Transformation (FFT) algorithm.

Now all that is left is a set of ordinary differential Eqs. 29, 30 and 31,
which can be solved with methods like finite-differences. We instead use the
Chebyshev method and discretize $y\in[0,1]$ to $M+1$ Chebyshev nodes, so
$y_{m}=[1+\cos(m\pi/M)]/2$ with $m=0,1,2,\dots M$. An advantage of using the
Chebyshev method is that the unevenly distributed Chebyshev nodes have a
higher resolution near the boundary $y=0$ and $y=1$, therefore resolving the
boundary layer structures. The differentiation operator $\frac{d}{dy}$ can be
approximated by a discrete operator $\mathbf{D}$ (Trefethen, 2000; Peyret,
2002), whose elements are

$\displaystyle D_{j,k}$
$\displaystyle=2\frac{c_{j}}{c_{k}}\frac{(-1)^{j+k}}{\cos(j\pi/M)-\cos(k\pi/M)}$
$\displaystyle 0\leq j,\,k\leq M,\,\,j\neq k,$ (36) $\displaystyle D_{j,j}$
$\displaystyle=-\frac{\cos(j\pi/M)}{1-\cos^{2}(j\pi/M)},$ $\displaystyle 1\leq
j\leq M-1,$ (37) $\displaystyle D_{0,0}$
$\displaystyle=-D_{M,M}=\frac{2M^{2}+1}{3}.$ (38)

Here $c_{0}=c_{M}=2$, $c_{j}=1$ for $1\leq j\leq M-1$.

The discrete operation on the LHS of Eq. 29 can therefore be written as

$\mathbf{A}=\mathbf{D}^{2}-\sigma_{k}\mathbf{I},$ (39)

where $\mathbf{I}$ is a $M+1$ by $M+1$ identity matrix.

Noticing that the Chebyshev nodes are
$\mathbf{y}=[1,y_{1},y_{2},\dots,y_{M-1},0]^{\top}$, we can write the discrete
solution to Eqs. 29, 30 and 31 at these locations as a column vector,

$\mathbf{U}=[\hat{g}_{1,k},\hat{u}_{k}(y_{1}),\hat{u}_{k}(y_{2}),\dots,\hat{u}_{k}(y_{M-1}),\hat{g}_{0,k}]^{\top}=[\hat{g}_{1,k},\tilde{\mathbf{U}},\hat{g}_{0,k}]^{\top}.$
(40)

The interior solution $\tilde{\mathbf{U}}$ therefore satisfies

$\tilde{\mathbf{A}}\tilde{\mathbf{U}}=\tilde{\mathbf{H}},$ (41)

where the $(M-1)\times(M-1)$ matrix $\tilde{\mathbf{A}}$ is the interior of
$\mathbf{A}$ (by removing its first and last rows and columns), and

$\displaystyle\tilde{\mathbf{H}}$
$\displaystyle=\begin{bmatrix}\hat{h}_{k}(y_{1})-\hat{g}_{1,k}A_{1,1}-\hat{g}_{0,k}A_{1,M}\\\
\hat{h}_{k}(y_{2})-\hat{g}_{1,k}A_{2,1}-\hat{g}_{0,k}A_{2,M}\\\ \vdots\\\
\hat{h}_{k}(y_{M-1})-\hat{g}_{1,k}A_{M-1,1}-\hat{g}_{0,k}A_{M-1,M}\end{bmatrix}.$
(42)

Equation 41 is invertible, and the operator $\tilde{\mathbf{A}}$ does not
change during time stepping while all the boundary information $g$ and forcing
$h$ are contained in $\tilde{\mathbf{H}}$. This allows us to compute the LU
decomposition of $\tilde{\mathbf{A}}$ at the beginning so
$\tilde{\mathbf{A}}\tilde{\mathbf{U}}=\tilde{\mathbf{H}}$ can be efficiently
inverted through backward and forward substitutions during time stepping. With
$\tilde{\mathbf{U}}$ and $\mathbf{U}$ solved, the Fourier coefficients
$\hat{u}_{k}$ can then be inserted into Eq. 28 and the solution $u(x,y)$ is
therefore obtained.

###  3.3 Helmholtz solver for Robin boundary conditions

Next we consider the Helmholtz solver for equations like Eq. 18, where
inhomogeneous Robin boundary conditions like Eq. 11 are applied. In a general
form, consider

$\displaystyle\nabla^{2}u-\sigma u$ $\displaystyle=h,$ (43) $\displaystyle
a_{0}(x)u(x,0)+b_{0}(x)$ $\displaystyle u_{y}(x,0)=g_{0}(x),$ (44)
$\displaystyle a_{1}(x)u(x,1)+b_{1}(x)$ $\displaystyle u_{y}(x,1)=g_{1}(x).$
(45)

The idea of solving these equations is to use the influence matrix method
(Peyret, 2002): A PDE with inhomogeneous Robin boundary conditions can be
converted into a series of PDEs with homogeneous Dirichlet boundary
conditions, which can be solved by the method detailed in Section 3.2.

We first separate the solution into several subproblems, so

$u(x,y)=\tilde{u}(x,y)+\sum_{l=1}^{2L}\xi_{l}\bar{u}_{l}(x,y),$ (46)

where $\xi_{l}$ are unknown coefficients to be determined later and

$\displaystyle\nabla^{2}\tilde{u}-\sigma\tilde{u}$ $\displaystyle=h,$ (47)
$\displaystyle\tilde{u}(x,0)=\tilde{u}(x,1)$ $\displaystyle=0,$ (48)
$\displaystyle\nabla^{2}\bar{u}_{l}-\sigma\bar{u}_{l}$ $\displaystyle=0,$ (49)
$\displaystyle\bar{u}_{l}|_{\eta_{m}}=\delta_{lm}\mbox{ for all }$
$\displaystyle\eta_{m}\in\partial\Omega.$ (50)

Here $\eta_{m}$ represents the $m$th node on the boundary, so there are $2L$
of them, and $\delta_{lm}$ is the Kronecker delta function. Now, the boundary
conditions Eqs. 44 and 45 indicate

$\left[a_{i}(x)\left(\sum_{l=1}^{2L}\xi_{l}\bar{u}_{l}\right)+b_{i}(x)\left(\frac{\partial\tilde{u}}{\partial
y}+\sum_{l=1}^{2L}\xi_{l}\frac{\partial\bar{u}_{l}}{\partial
y}\right)-g_{i}(x)\right]\biggr{\rvert}_{\eta_{m}}=0\mbox{ for
}\eta_{m}\in\partial\Omega_{i}.$ (51)

Here $i\in\\{0,1\\}$ indicates the bottom or top boundary. As Eq. 51 holds on
all the boundary nodes $\eta_{m}$, it provides $2L$ equations for $2L$
unknowns $\xi_{l}$, and such a linear system is invertible. With all $\xi_{l}$
solved, the solution of Eqs. 43, 44 and 45 can be recovered as
$u=\tilde{u}+\sum_{l=1}^{2L}\xi_{l}\bar{u}_{l}$.

In fact, $\xi_{m}$ is exactly the Dirichlet data for $u$ at boundary node
$\eta_{m}$, therefore the solution of Eqs. 43, 44 and 45 is the same as the
solution of

$\displaystyle\nabla^{2}u-\sigma u$ $\displaystyle=h,$ (52) $\displaystyle
u|_{\eta_{m}}$ $\displaystyle=\xi_{m}\mbox{ for }\eta_{m}\in\partial\Omega.$
(53)

This method can directly solve the heat Eq. 18 with boundary conditions in
Eqs. 8 and 12, by assigning $a_{0}=1$, $b_{0}=0$, $g_{0}=1$ at the bottom and
$a_{1}=1-a$, $b_{1}=a$, $g_{1}=0$ at the top. Usually, the solutions of Eqs.
49 and 50 are obtained and saved at the beginning. At each time step, Eqs. 47
and 48 are solved and the location of moving plates will determine $a_{i}(x)$,
$b_{i}(x)$, and $g_{i}(x)$, so Eq. 51 can be inverted to provide $\xi_{m}$
which can be used as boundary data in Eqs. 52 and 53, and the solution
$\theta_{n}$ can therefore be determined.

###  3.4 Flow solver with Robin boundary conditions

With the heat Eq. 18 solved, the only missing part is the following flow
problem at time step $n$,

$\displaystyle\nabla^{2}\omega-\sigma\omega$ $\displaystyle=f,$ (54)
$\displaystyle-\nabla^{2}\psi$ $\displaystyle=\omega,$ (55)
$\displaystyle\psi=\psi_{y}$ $\displaystyle=0\quad\quad\mbox{at }y=0,$ (56)
$\displaystyle\psi=0,\quad a\,\psi_{y}+(1-a)\,\psi_{yy}$
$\displaystyle=g\quad\quad\mbox{at }y=1.$ (57)

Here we have dropped all the subscripts. We can also solve these equations
with the influence matrix method. Now we want to convert the Neumann boundary
condition in Eq. 56 and Robin boundary condition in Eq. 57 to a Dirichlet
boundary condition for vorticity $\omega$. We decompose $\omega$ and $\psi$
as,

$\displaystyle\omega(x,y)=\tilde{\omega}(x,y)+\sum_{l=1}^{2L}\xi_{l}\bar{\omega}_{l}(x,y),$
(58)
$\displaystyle\psi(x,y)=\tilde{\psi}(x,y)+\sum_{l=1}^{2L}\xi_{l}\bar{\psi}_{l}(x,y).$
(59)

The associated subproblems are

$\displaystyle\nabla^{2}\tilde{\omega}-\sigma\tilde{\omega}$
$\displaystyle=f,\,\,-\nabla^{2}\tilde{\psi}=\tilde{\omega}$ (60)
$\displaystyle\tilde{\omega}(x,0)=\tilde{\omega}(x,1)$
$\displaystyle=\tilde{\psi}(x,0)=\tilde{\psi}(x,1)=0$ (61)
$\displaystyle\nabla^{2}\bar{\omega}_{l}-\sigma\bar{\omega}_{l}$
$\displaystyle=0,\,\,-\nabla^{2}\bar{\psi}_{l}=\bar{\omega}_{l}$ (62)
$\displaystyle\bar{\psi}_{l}(x,1)=0,\,\,\bar{\omega}_{l}|_{\eta_{m}}$
$\displaystyle=\delta_{lm}\mbox{ for all }\eta_{m}\in\partial\Omega.$ (63)

The Neumann boundary condition in Eq. 56 and the Robin boundary condition in
Eq. 57 can now be enforced as

$\displaystyle\left(\frac{\partial\tilde{\psi}}{\partial
y}+\sum_{l=1}^{2L}\xi_{l}\frac{\partial\bar{\psi}_{l}}{\partial
y}\right)\biggr{\rvert}_{\eta_{m}}=0$ $\displaystyle\mbox{ for
}\eta_{m}\in\partial\Omega_{0},$ (64)
$\displaystyle\left[a\left(\frac{\partial\tilde{\psi}}{\partial
y}+\sum_{l=1}^{2L}\xi_{l}\frac{\partial\bar{\psi}_{l}}{\partial
y}\right)+(1-a)\left(\frac{\partial^{2}\tilde{\psi}}{\partial
y^{2}}+\sum_{l=1}^{2L}\xi_{l}\frac{\partial^{2}\bar{\psi}_{l}}{\partial
y^{2}}\right)-g\right]\biggr{\rvert}_{\eta_{m}}=0$ $\displaystyle\mbox{ for
}\eta_{m}\in\partial\Omega_{1}.$ (65)

Equations 64 and 65 are again $2L$ equations for $2L$ unknowns $\xi_{l}$ and
the linear system is invertible. Equations 58 and 59 can then be summed and
the solutions are obtained. In fact, $\xi_{l}$ is exactly the boundary data
for $\omega$ at boundary node $\eta_{l}$, so we can instead solve

$\displaystyle\nabla^{2}\omega-\sigma\omega$ $\displaystyle=f,$ (66)
$\displaystyle-\nabla^{2}\psi$ $\displaystyle=\omega,$ (67)
$\displaystyle\psi|_{\eta_{m}}=0,\,\,\omega|_{\eta_{m}}=\xi_{m}$
$\displaystyle\mbox{ for }\eta_{m}\in\partial\Omega.$ (68)

###  3.5 Smooth boundary conditions

In principle, the introduced fluid and heat solver is able to manage the
inhomogeneous Robin boundary condition Eq. 11 at $y=1$. However this boundary
condition is not smooth, therefore limiting the accuracy of a numerical method
with finite resolution. To overcome this, we aim to construct a smooth
characteristic function $\hat{\mathbb{1}}_{x\in P_{i}}$ so it is compactly
supported and sufficiently smooth.

We first construct a smooth step function in 1D, whose derivative
$\phi_{l,m}(r)$ is in the family of Wendland functions that are shaped like a
Gaussian (Chernih et al., 2014):

$\phi_{l,m}(r)=\left\\{\begin{array}[]{cl}\frac{1}{\Gamma(m)2^{m-1}}\int_{r}^{1}s(1-s)^{l}(s^{2}-r^{2})^{m-1}ds&\mbox{for
}0\leq r\leq 1,\\\ 0&\mbox{for }r>1.\end{array}\right.$ (69)

The integer $m$ controls the smoothness of the Wendland function, and
$l=\left\lfloor m+n/2\right\rfloor+1$ for spatial dimension $n$. It can be
shown that $\phi_{l,m}\in C^{2m}(\mathbb{R}^{+})$ and it is compactly
supported. Next, we take $m=1$, $l=2$ and construct a smooth step function
$W_{\epsilon}(x)$ that transitions from 0 to 1 on $[-\epsilon,\epsilon]$,

$W_{\epsilon}(x)=\frac{\int_{-\infty}^{x}\phi_{2,1}(|s|/\epsilon)ds}{\int_{-\infty}^{\infty}\phi_{2,1}(|s|/\epsilon)ds}.$
(70)

This function is plotted in Fig. 3(a) with various $\epsilon$, and it is easy
to verify that $W_{\epsilon}(x)=0$ when $x<-\epsilon$ and $W_{\epsilon}(x)=1$
when $x>\epsilon$. As the transition length $2\epsilon$ becomes smaller,
$W_{\epsilon}$ becomes sharper, approximating a step function. Moreover,
$W_{\epsilon}\in C^{3}(\mathbb{R})$ when $\epsilon>0$ due to our choice of
$m=1$, achieving our goal of constructing a smooth step function.

Figure 3: Smooth step and indicator functions. (a) Smooth step function
$W_{\epsilon}$ that has a transition interval of $[-\epsilon,\epsilon]$. Four
values of $\epsilon=0.01,0.1,0.2,0.4$ are plotted. (b) Smooth indicator
function $\hat{a}$ for locating the region of solid plates. The parameters
plotted are $x_{p}^{(1)}=1$, $x_{p}^{(2)}=3$, $d=1$, and $\epsilon=0.05$.

We can further evaluate $W_{\epsilon}(x)$ as

$W_{\epsilon}(x)=\begin{cases}0\quad\mbox{if }x<-\epsilon,\\\
-\frac{3}{4}\left(\frac{x}{\epsilon}\right)^{5}+\frac{5}{2}\left(\frac{x}{\epsilon}\right)^{4}\mbox{sgn}(x)-\frac{5}{2}\left(\frac{x}{\epsilon}\right)^{3}+\frac{5}{4}\frac{x}{\epsilon}+\frac{1}{2}\quad\mbox{if
}x\in[-\epsilon,\epsilon],\\\ 1\quad\mbox{if }x>\epsilon.\\\ \end{cases}$

We next construct a smooth characteristic function $\hat{\mathbb{1}}_{x\in P}$
for a plate $P$ centered at $x_{p}$ with length $d$,

$\hat{\mathbb{1}}_{x\in
P_{i}}=W_{\epsilon}\left(x-(x_{p}-\frac{d}{2})\right)W_{\epsilon}\left((x_{p}+\frac{d}{2})-x\right).$
(71)

In our numerical examples, we take $\epsilon=0.05d$ to ensure smoothness.
Finally, we can write the Robin boundary conditions in the following form.

$\begin{cases}(1-\hat{a})\,\theta+\hat{a}\,\theta_{y}=0\\\
\hat{a}\,u+(1-\hat{a})\,u_{y}=\hat{g}\quad\quad\mbox{at }y=1,\\\
v=0\end{cases}$ (72)

where

$\hat{a}(x)=\sum_{i}\hat{\mathbb{1}}_{x\in
P_{i}},\,\,\hat{g}(x)=\sum_{i}u_{p}^{(i)}\hat{\mathbb{1}}_{x\in P_{i}}.$ (73)

Noticing that $W_{\epsilon}(0)=0.5$, it can be easily verified that
$\hat{a}(x)\in[0,1]$ as long as $|x_{p}^{(i)}-x_{p}^{(j)}|\geq d$ for $i\neq
j$. The function $\hat{a}(x)$ for two plates centered at $x_{p}^{(1)}=1$,
$x_{p}^{(2)}=3$ with plate length $d=1$ is shown in Fig. 3(b).

###  3.6 Dynamics of the moving plate

There are two types of forcing on each plate. One is the fluid force due to
shear stress, the other is the interaction force when two plates make contact.

For the fluid force, we simply integrate the shear stress by replacing the
characteristic function with its smooth version,

$f^{(i)}=-\text{Pr}{}\int_{0}^{\Gamma}u_{y}(x,1,t)\hat{\mathbb{1}}_{x\in
P_{i}}d\,x=-\text{Pr}{}\Gamma\langle u_{y}\hat{\mathbb{1}}_{x\in
P_{i}}\rangle.$ (74)

As we are using an equally spaced periodic grid in $x$, the integration can be
replaced with a numerical average of all grid values of the integrand so
$\int_{0}^{\Gamma}f(x)\,dx\approx\Gamma\langle
f\rangle=(\Gamma/L)\sum_{k=1}^{L}f_{k}$, which is spectrally accurate
(Trefethen, 2000).

When solid plates make contact, a contact force between them keeps the plates
separated. Numerically, we approximate this force such that it is short-ranged
and stays smooth until two plates are too close to each other. For plate
number $i$, there are two forces for contact from the left neighbor and from
the right neighbor,

$\displaystyle f_{l}^{(i)}$
$\displaystyle=f_{\max}W_{\delta}\left(d-|x_{p}^{(i)}-x_{p}^{(i-1)}|\right),$
(75) $\displaystyle f_{r}^{(i)}$
$\displaystyle=-f_{\max}W_{\delta}\left(d-|x_{p}^{(i+1)}-x_{p}^{(i)}|\right).$
(76)

Here the parameter $\delta$ models an “interaction length”, and $f_{\max}$ is
the maximum interacting force between two plates. At each simulation, $\delta$
and $f_{\max}$ are chosen according to the spatial and temporal resolution, so
the ODE and PDE solvers can sufficiently resolve the plate motion and the
associated boundary conditions. We note that this choice of interacting force
conserves both the total momentum and kinetic energy of plates when fluid
force is not present.

###  3.7 Summary of numerical methods

There are two stages during our numerical simulation and below we list some of
the key steps during each stage.

Preprocessing stage:

  1. 1.

The inverse $\tilde{\mathbf{A}}^{-1}$ (or the LU decomposition of
$\tilde{\mathbf{A}}$) in Eq. 41 is prepared for operators in Eqs. 17, 18 and
19, by taking $\sigma=0,\sigma_{1},\sigma_{2}$ in Eq. 39;

  2. 2.

Subproblems Eqs. 49, 50, 62 and 63 are solved, and the solutions
$\bar{\theta}_{l}$, $\bar{\omega}_{l}$, and $\bar{\psi}_{l}$ for
$l=1,2,\dots,2L$ are saved.

At the $n$th step of the time-stepping stage:

  1. 1.

The fluid and interaction forces on each plate are computed according to Eqs.
74, 75 and 76;

  2. 2.

The location and velocity of each plate are evolved by Eqs. 23 and 24, and the
indicator function $\hat{a}$ is prepared by Eq. 73;

  3. 3.

The forcing term $h_{n}$ in Eq. 18 is prepared according to Eq. 22, and
$\tilde{\theta}_{l}$ is solved from Eqs. 47 and 48;

  4. 4.

Equation 51 is inverted so the Dirichlet boundary data for $\theta_{n}$ is
known, Eqs. 52 and 53 are then solved for $\theta_{n}$;

  5. 5.

The temperature $\theta_{n}$ is used to prepare $f_{n}$ according to Eq. 21,
so $\tilde{\omega}_{l}$, $\tilde{\psi}_{l}$ can be solved from Eqs. 60 and 61;

  6. 6.

Equations 64 and 65 are inverted so the Dirichlet data for $\omega_{n}$ is
known, Eqs. 67 and 68 are finally solved to provide $\omega_{n}$ and
$\psi_{n}$.

##  4 Results

###  4.1 One plate dynamics

Figure 4:  Motion of a small plate ($\mbox{{Cr}}{}=0.1$) is random and
bidirectional. (a) A snapshot of flow and temperature fields beneath a plate.
The small plate is trapped at a cool converging center. (b) Vertically
averaged temperature $\bar{\theta}$ and vertical velocity $\bar{v}$ at the
same moment of (a). The shaded region indicates the location of the plate. At
the converging center, the averaged temperature is low and the flow moves
downward. (c)-(d) The displacement $x_{p}$ and velocity $u_{p}$ of the plate
show behavior of a random walk with jumps.

In this section, we review the dynamics of a single plate motion. To simplify
our study, the Rayleigh number is fixed at $\text{Ra}=10^{6}$, the Prandtl
number is $\text{Pr}=7.9$, and the aspect ratio is $\Gamma=4$. These
parameters are similar to previous numerical study (Huang, 2024). For the
plate, we set the mass as $m=4d$, so plates with various length $d$ have the
same density. For the numerical solver, there are $512$ Fourier modes in the
$x$ direction and $129$ Chebyshev nodes in the $y$ direction, and the time
step size is $\Delta t=10^{-6}$. These parameters yield accurate, stable, and
resolved numerical solutions.

To measure the size of the plate, we define the covering ratio
$\mbox{{Cr}}{}=d/\Gamma$. Depending on the size of the plate, or Cr, the
dynamics of the plate motion can be very different. Figure 4 and Supplemental
Movie 1 show the dynamics of a small plate with $\mbox{{Cr}}{}=0.1$, and its
motion is a continuous random walk shown in Fig. 4(c). However in Fig. 5 and
Supplemental Movie 2, a larger plate with $\mbox{{Cr}}=0.6$ shows completely
different dynamics: it translates unidirectionally as shown in Fig. 5(c).

In Huang (2024), the two dynamics are analyzed in detail, and we summarize the
key interplay between the plate and the fluid below.

In Fig. 4, the small plate tends to be attracted by the converging center of
the fluid – the location where the fluid sinks. This converging center can be
clearly seen in Fig. 4(a)-(b), located at the minimum of both the vertically
averaged temperature $\bar{\theta}(x,t)=\int_{0}^{1}\theta(x,y,t)\,dy$, and
the averaged vertical flow speed $\bar{v}(x,t)=\int_{0}^{1}v(x,y,t)\,dy$. This
means the plate velocity $u_{p}=\dot{x}_{p}$ in Fig. 4(d) matches the
translational velocity of the flow converging center, which has a 0 mean but
is subject to noise due to random fluid forcing. In this case, the plate
motion is passive and completely driven by the flow structure, and the
converging center of the surface flow serves as a stable equilibrium location
for the plate.

Figure 5:  Motion of a large plate ($\mbox{{Cr}}{}=0.6$) is unidirectional.
(a)-(b) Flow and temperature fields beneath the plate. (c)-(d) The
displacement $x_{p}$ and velocity $u_{p}$ of the moving plate, which shows a
unidirectional motion with nonzero mean velocity.

For a larger plate with $\mbox{{Cr}}{}=0.6$, Fig. 5 and Supplementary Movie 2
show that the plate motion becomes unidirectional. Increasing plate size
clearly changes the flow and temperature distribution in the fluid, as the
bulk fluid temperature in Fig. 5(a)-(b) is visibly higher than that in Fig.
4(a)-(b). This is a clear sign of the thermal blanket effect, as the bigger
plate shields the heat from escaping and the effective cooling area at $y=1$
is smaller. In this case, the plate is no longer passive, but creates a
thermal blanket that warms the fluid beneath it. Unlike the situation of small
plates, a large plate sitting on top of a converging center cannot be stable
in the long term, as eventually the temperature beneath the plate will become
high enough to turn this converging center into a divergent one. Shown in Fig.
5(b), the average temperature $\bar{\theta}$ is indeed higher below the plate,
and the plate sits between the converging and diverging centers. This causes
the unidirectional motion of the plate, and as the plate keeps affecting the
temperature distribution beneath it, the temperature and flow fields move with
the plate as shown in Supplementary Movie 2. The plate displacement $x_{p}$
and velocity $u_{p}$ are shown in Fig. 5(c)-(d), where $u_{p}$ has a nonzero
mean and is subject to random forcing from the fluid.

Figure 6: Plate displacement and velocity for different covering ratios Cr.
(a) Sample trajectories of the plate location, where small plates are more
affected by noise and large plates have more persistent unidirectional motion.
(b) Total travel of the plate reveals its speed, a maximum speed of travel can
be seen at around $\mbox{{Cr}}{}=0.6$. (c) Average travel speed has a maximum
at $\mbox{{Cr}}{}=0.56$, and unidirectional motions start to appear for plates
larger than $\mbox{{Cr}}{}=0.33$.

The motions of plates with various Cr are shown in Fig. 6. The displacement in
Fig. 6(a) clearly shows that the small plate has a random motion whose net
displacement grows slowly in time. As Cr increases, the plate starts to have
more persistent unidirectional motions, however the random fluid forcing can
easily reverse the travel direction of the plate, leading to reversals of
direction in Fig. 6(a). Further increasing Cr, the unidirectional motion
becomes more persistent, however the velocity (slope of $x_{p}$) decreases.
Defining the total travel of a plate as
$d_{p}(t)=\int_{0}^{t}|u_{p}(t^{\prime})|\,dt^{\prime}$, Fig. 6(b) shows a
peak of the plate traveling speed at around $\mbox{{Cr}}{}=0.6$. To further
verify this, we define the average plate speed as
$v_{p}=\lim_{t\to\infty}d_{p}(t)/t$ in Fig. 6(c), and a maximum indeed appears
at $\mbox{{Cr}}{}=0.6$.

Figure 7:  Time series (lower panels) and histogram (upper panels) of the
plate velocity $u_{p}$ at various Cr. The plate velocity is normalized by its
average travel speed $v_{p}$, so $u_{p}\approx\pm v_{p}$ suggests a
unidirectional translation. (a)-(e) have the covering ratio of $0.0625$,
$0.3125$, $0.375$, $0.4375$, $0.875$. The plate motion has a transition from
the passive state (a)-(b) to the translating state (d)-(e), and the
translation is also more persistent for large Cr in (e).

As the thermal blanket effect strengthens, there is an apparent transition of
the plate dynamics. Figure 7 shows the time series (lower panels) and
histogram (upper panels) of $u_{p}$ at various Cr. For small Cr, the histogram
of plate velocity resembles a Gaussian distribution, whose zero mean suggests
that the net plate displacement would be small. Increasing Cr beyond $0.3$,
the plate dynamics start to transition as Fig. 7$(c)$ shows the variance of
$u_{p}$ increases. At $\mbox{{Cr}}=0.4375$ [Fig. 7(d)], the two translational
states with $u_{p}=\pm v_{p}$ emerge, where $u_{p}$ switches between the two
directions due to the noise of fluid forcing. At even higher Cr [Fig. 7(e)],
the unidirectional motion is persistent and the reversal becomes rare. The
observation here matches the stochastic theory developed in Huang (2024),
which consists of a simple model that recovers the mechanical and thermal
interplay between the plate and the fluid. This stochastic model predicts that
there is only a passive state (no net plate motion) for
$\mbox{{Cr}}<\mbox{{Cr}}^{*}$ where the critical covering ratio is
$\mbox{{Cr}}^{*}=1/3$ for $\Gamma=4$, and the translational states can only
exist for plates with $\mbox{{Cr}}>\mbox{{Cr}}^{*}$, indeed matching Fig. 7.

Figure 8:  Nusselt and Reynolds numbers for the convecting flow. (a) The
Nusselt number is a measure of the vertical heat passing through the flow. (b)
Reynolds number as a measure of flow speed. The free data are for the plate
freely moving with the flow, and the immobile data are for the plate that is
fixed.

Finally, we investigate how the bulk properties of the flow respond to the
moving plate. In Fig. 8, we show the Nusselt number
$\text{Nu}=-[\Gamma(t_{2}-t_{1})]^{-1}\int_{t_{1}}^{t_{2}}dt\int_{0}^{\Gamma}\theta_{y}(x,0,t)\,dx$,
and Reynolds number
$\text{Re}=(t_{2}-t_{1})^{-1}\int_{t_{1}}^{t_{2}}\left(\max_{(x,y)}|\bm{u}(x,y,t)|\right)\,dt$,
where $t_{2}-t_{1}$ is the interval for long-time average. The two groups of
measurements are for a plate that is free to move by the flow (free), and for
a plate that is fixed at a certain location (immobile). By setting the plate
free, the Nusselt number decreases slightly, and this is possibly due to the
moving plate effectively shielding a larger portion of the free surface. With
a stronger shielding effect, the vertical heat transfer measured by the
Nusselt number decreases. The flow speed, on the other hand, is seen to
decrease once we let the plate move freely. This can either be caused by the
increased shielding of the free surface, or be a consequence of the flow
passing its kinetic energy to the moving plate. We also note at the critical
$\mbox{{Cr}}^{*}$, the flow speed reaches a maximum for the immobile plate
shown in Fig. 8(b), and the Nusselt number approaches its maximum at around
$\mbox{{Cr}}=0.6$ where the plate translates the fastest as shown in Fig.
6(c).

###  4.2 Two-plate interactions

Figure 9:  Dynamics of two small plates ($\mbox{{Cr}}{}_{p}=0.1$ each) forming
a supercontinent of $\mbox{{Cr}}=0.2$. (a) Flow and temperature distribution
beneath the supercontinent. The surface flow is converging and the formation
of the supercontinent is stable. (b) Vertically averaged temperature
$\bar{\theta}$ and vertical velocity $\bar{v}$ at the same moment of (a), with
the region of the two plates shaded. (c) The displacement of plate
$x_{p}^{(1)}$ and $x_{p}^{(2)}$. (d) The normalized plate distance
$d_{12}=[x_{p}^{(2)}-x_{p}^{(1)}]/\Gamma$ indicates the two plates tend to
stay in contact. The white region (plates separated) and gray region (plates
in contact) are separated by $\mbox{{Cr}}_{p}$ and $1-\mbox{{Cr}}_{p}$.

Adding multiple plates to the convective surface brings interactions between
plates and leads to more diverse dynamics. In our numerical simulation, adding
a second plate can easily be achieved through the indicator function method
outlined in Sections 2 and 3. In the following numerical experiments, we set
$\text{Ra}=10^{6}$, $\Pr=7.9$, $\Gamma=4$, $m=4d$ as we did for the single
plate case. We additionally set the maximum interaction force
$f_{\max{}}=10^{6}$ and an interaction range $\delta=\epsilon$ that matches
the size of smoothing region of the indicator function in Section 3.5. These
two parameters define the force of interaction between the two plates through
Eqs. 75 and 76, and such interaction conserves both the kinetic energy and
momentum of the plates.

The dynamics of a pair of small plates (Fig. 9) and a pair of large plates
(Fig. 10) are quite different. In Fig. 9 and Supplementary Movie 3, two small
plates with individual covering ratios $\mbox{{Cr}}_{p}=0.1$ are released on
the convective surface. The two plates tend to stay together, generating a
“supercontinent” as shown in Fig. 9(a). Further analyzing the flow temperature
and surface flow rate in Fig. 9(b), we see they are in fact attracted by a
converging center of the surface flow, and the surface flow pushes them into
each other. The trajectories $x_{p}^{(1)}$ and $x_{p}^{(2)}$ of the plates are
shown in Fig. 9(c), and the normalized plate distance
$d_{12}=[x_{p}^{(2)}-x_{p}^{(1)}]/\Gamma$ is plotted in Fig. 9(d). We clearly
see that the two plates prefer to stay in contact, as the normalized distance
stays near $\mbox{{Cr}}_{p}$ or $1-\mbox{{Cr}}_{p}$ in Fig. 9(d).

The combined covering ratio of these two plates is
$\mbox{{Cr}}=2\mbox{{Cr}}_{p}=0.2$, less than the critical covering ratio of
$\mbox{{Cr}}^{*}=1/3$ we identified earlier. Therefore the thermal blanket
effect generated by this supercontinent is not strong enough to heat up the
fluid beneath, and the surface flow stays converging and pushing the two
plates together. Thus, a supercontinent with combined
$\mbox{{Cr}}<\mbox{{Cr}}^{*}$ is stable in its formation and exhibits a
passive motion.

Figure 10 and Supplementary Movie 4 show the dynamics of two plates with
$\mbox{{Cr}}=2\mbox{{Cr}}_{p}=0.6$. In this case, the fluid beneath the
supercontinent is warmed up due to the thermal blanket effect and generates an
upwelling flow. The resulting diverging surface flow separates the two plates,
leading to an unstable supercontinent formation. Figure 10(c)-(d) show the
plate trajectories and the normalized plate distance, and the two plates are
seen to stay separated in Fig. 10(d) as their normalized distance is in
between $\mbox{{Cr}}_{p}$ and $1-\mbox{{Cr}}_{p}$ and the contact sate is only
transient.

Figure 10: Two large plates ($\mbox{{Cr}}{}_{p}=0.3$ each) cannot form a
stable supercontinent of $\mbox{{Cr}}=2\mbox{{Cr}}{}_{p}=0.6$. (a) The warm
upwelling fluid creates a diverging surface flow beneath the plates. This
diverging surface flow pulls the two plates apart, rendering the
supercontinent formation unstable. (b) Vertically averaged temperature and
vertical velocity of the fluid beneath the plates shown in (a). (c) Plate
displacement $x_{p}^{(1)}$ and $x_{p}^{(2)}$. (d) The normalized plate
distance $d_{12}$ stays in the white region where the two plates are
separated.

In Fig. 10, although the covering ratio for each plate
$\mbox{{Cr}}_{p}<\mbox{{Cr}}^{*}$, their combined
$\mbox{{Cr}}=2\mbox{{Cr}}_{p}>\mbox{{Cr}}^{*}$. The supercontinent, once
formed, will become unstable as the warm fluid beneath creates a diverging
surface flow that pulls the two plates apart.

We now see the role of $\mbox{{Cr}}^{*}$ in the formation of supercontinents:
It is only possible to have a stable supercontinent if its Cr is less than
$\mbox{{Cr}}^{*}$.

Figure 11:  Contact and motion of the plates depend on the covering ratio. (a)
The normalized contact time $\tau$ decreases sharply when Cr increases above
$\mbox{{Cr}}^{*}=1/3$. (b) The plate center of mass velocity $v_{\mbox{COM}}$
increases when $\mbox{{Cr}}>\mbox{{Cr}}^{*}$, indicating the plates are no
longer passive to the flow.

Figure 11 shows how the plate dynamics depend on the combined Cr. In Fig.
11(a), we define a normalized contact time $\tau=t_{c}/T$, where $t_{c}$ is
the amount of time that the two plates are in contact and $T$ is the total
simulation time. A sharp decrease of $\tau$ is seen near
$\mbox{{Cr}}^{*}=1/3$, beyond which a persistent formation of supercontinent
is unlikely. This is consistent with our analysis earlier, as a supercontinent
with $\mbox{{Cr}}>\mbox{{Cr}}^{*}$ would induce warm upwelling flows that
disintegrate the supercontinent. The center of mass (CoM) velocity of plates
also picks up when $\mbox{{Cr}}>\mbox{{Cr}}^{*}$ [Fig. 11(b)], suggesting the
two large plates are no longer passive but instead translate like we have seen
in Fig. 6(c).

###  4.3 Multiple plates

Further increasing the number of plates, the formation of supercontinents
exhibits complex and intriguing dynamics. In Fig. 12 and Supplementary Movie
4, 8 plates with $\mbox{{Cr}}_{p}=0.057$ (total $\mbox{{Cr}}=0.457$) are
released on top of the convective fluid, where $\text{Ra}=10^{6}$, $\Pr=7.9$,
$\Gamma=4$, $m=4d$. The plate maximum interaction force $f_{\max{}}=10^{6}$
and interaction range $\delta=\epsilon$ are the same as before.

Figure 12(a) shows a moment when the 8 plates form 2 supercontinents, with
each supercontinent covering $4\mbox{{Cr}}_{p}=0.23$ of the free surface. Each
supercontinent thus has a covering ratio less than $\mbox{{Cr}}^{*}=1/3$ and
is thereby stable by our earlier argument. Indeed, the stable formation of a
supercontinent of 4 plates can be seen on the plate trajectory
$(x_{p}^{(1)},x_{p}^{(2)},\cdots,x_{p}^{(8)})$ shown in Fig. 12(b), whose
zoomed-in view between $t\in(0.35,0.36)$ is shown in Fig. 12(c). Of course,
the formation of a supercontinent with a different number of plates is
possible, and our theory predicts that they are stable if the number of plates
$I<\mbox{{Cr}}^{*}/\mbox{{Cr}}_{p}\approx 6$.

Figure 12:  Dynamics of 8 plates ($\mbox{{Cr}}{}_{p}=0.057$,
$\mbox{{Cr}}=8\mbox{{Cr}}{}_{p}=0.457$) floating on top of the convecting
fluid. (a) A snapshot of 8 plates and the convective fluid beneath, the center
location of the plates are $(x_{p}^{(1)},x_{p}^{(2)},\cdots,x_{p}^{(8)})$. (b)
Trajectories of $(x_{p}^{(1)},x_{p}^{(2)},\cdots,x_{p}^{(8)})$, and plates can
be seen forming supercontinents over time. (c) Zoomed-in view of the
trajectories in (b) in the time window of $t\in(0.35,0.36)$.

To further analyze the formation of supercontinents, we define the formation
number $I(t)$ as the number of plates forming the largest supercontinent at
any given time $t$. A schematic of such a formation number is shown in Fig.
13(a), where 5 continents are formed and the largest supercontinent has $I=4$.
The formation number can be tracked over time, and Fig. 13(b) shows the
formation number of the simulation presented in Fig. 12, with the zoomed-in
view of $I(t)$ in the window of Fig. 12(c) shown in Fig. 13(c). We note that
the formation number can actually take on all integer values between 1 and 8,
although many of the formation numbers are transient (such as $I=1$ and
$I=8$). The most common and persistent formation number we can visually see in
Fig. 13(c) is $I=4$, and this is confirmed in the histogram of $I(t)$ shown in
Fig. 13(d).

In Fig. 13(d), each bin corresponds to the total number of appearances of
supercontinents with size $I$ (shown on top of each bin), and the horizontal
axis shows the corresponding covering ratio. We indeed see that $I=4$ plates
is the most frequent size of a supercontinent. For small supercontinents, they
tend to merge and form a larger but stable supercontinent. For large
supercontinents, the random fluid forcing can trigger a disintegration and
break it into many smaller ones. Especially for $I\geq 6$, our theory predicts
such formations are not stable, hence supercontinents with $I\geq 6$ are
rarely formed and (once formed) are unstable.

Figure 13:  Statistics of the formation number $I(t)$ that is the maximum
number of plates in a supercontinent at time $t$. (a) Schematic of $I=4$. (b)
Time series of $I(t)$ shows the possibility of forming supercontinents with
various sizes. (c) Zoomed-in view of $I(t)$ in (b) during $t\in(0.35,0.36)$.
(d) Histogram of $I(t)$ indicates $I=4$ is the most common supercontinent
formation, while small and large supercontinents are unlikely to form. The
histogram is plotted against the size of supercontinent
$\mbox{{Cr}}=I\mbox{{Cr}}_{p}$ and the formation number $I$ is labelled on top
of each bin.

The results shown in Fig. 12 and Fig. 13 are common for simulations with
different number ($N_{p}$) and size ($\mbox{{Cr}}_{p}$) of plates. And we
identify $I\mbox{{Cr}}_{p}<\mbox{{Cr}}^{*}$ as a clue for predicting the
supercontinental size $I$.

##  5 Discussion

Through our numerical investigations, we clearly see the covering ratio as the
main factor affecting the thermal blanket effect, which determines the plate
dynamics. This is especially apparent during the formation of supercontinents–
the continent covering ratio cannot exceed $\mbox{{Cr}}^{*}$, as the strong
thermal blanket effect will induce a diverging surface flow that pulls the
formed supercontinent apart.

As our current study is inspired by laboratory experiments, we certainly look
forward to future experimental investigations of the interaction between
multiple plates. Besides the geometry presented in Figs. 1 and 2, the broader
investigation of fluid-structure interactions in thermal convection also
includes adding fixed obstacles to the convective domain (Bao et al., 2015; Li
et al., 2024), modifying the convective boundary conditions (Zhang et al.,
2020; Huang & Zhang, 2022), and allowing for moving boundaries that are driven
by the fluid forcing (Mercier et al., 2014; Wang & Zhang, 2023) or phase
change (Huang et al., 2021; Huang & Moore, 2022; Weady et al., 2022). Many of
these works are certainly geophysically inspired, addressing long-standing
mysteries such as the super rotation of Earth’s inner core (Wang & Zhang,
2023; Yang & Song, 2023).

In order to more closely simulate the geophysical plate tectonics, there are
several future directions to improve our model. The first one is to
investigate the plate tectonics in a 3-dimensional geometry, including the
mantle-like spherical shell and the rectangular cuboidal fluid domain that is
periodic in two horizontal directions. Pioneer works including Lowman & Gable
(1999); Mao et al. (2019); Mao (2021) have investigated the fluid-structure
interaction and plate interaction in these settings. They also used the
geophysical parameters of the mantle, most notably a Prandtl number around
$10^{23}$ (Meyers et al., 1987). Incorporating the aforementioned geometry and
parameters of the mantle is one future direction for our study.

The interaction between continental plates is also more complicated in
geophysical plate tectonics, as the converging continental plates deform the
contact region and form the tallest mountains of Earth. It is still an ongoing
quest to understand the consequences of converging and diverging continents,
with many recent works focusing on the deforming contact region and addressing
its influence on the mantle flow beneath (Rozel et al., 2017; Whitehead,
2023). Inspired by the geophysics, one modification to our current model could
be changing the way neighboring plates interact. This can be achieved through
changing the coefficient of restitution to 0 during each collision, which
better captures the collision between continental plates during a long time
scale. Adding an attractive force between plates can also reveal interesting
dynamics, as the diverging surface flow has to be strong enough to pull the
plates apart. These modifications are awaiting implementation, and they may
reveal new insights into the physics and dynamics of the formation of
supercontinents.

In our present work, we have fixed the fluid aspect ratio $\Gamma=4$ in all
the numerical examples. We note that $\Gamma$ can significantly affect the
formation of supercontinents. In Huang (2024), the critical covering ratio
$\mbox{{Cr}}^{*}$ is shown to depend on the aspect ratio as
$\mbox{{Cr}}^{*}\sim\Gamma^{-2/3}$ when $\Gamma$ is large. This indicates that
the maximum length of a stable supercontinent is
$\Gamma\mbox{{Cr}}^{*}\sim\Gamma^{1/3}$, so its size increases slowly with
increasing $\Gamma$.

To conclude, we consider a toy model to predict the size of the largest
continental plate. The aspect ratio of Earth’s mantle is roughly
$\Gamma_{E}=10$, so the critical covering ratio there is
$\mbox{{Cr}}^{*}_{E}\approx(\Gamma_{E}/\Gamma)^{-2/3}\mbox{{Cr}}^{*}\approx
0.2$ where we have taken $\Gamma=4$ and $\mbox{{Cr}}^{*}=1/3$ from our study.
This gives an estimated dimension of the largest stable continental plate of
Earth, $L\approx 2\pi R_{E}\mbox{{Cr}}^{*}_{E}$, where $R_{E}=6400$ km is the
radius of Earth. So the plate area is approximately $L^{2}\approx 7\times
10^{7}$ km2, slightly underestimating the largest continental plate, the North
American Plate, whose area is $7.59\times 10^{7}$ km2.

## Supplementary Materials

Supplementary movies are available at
https://math.nyu.edu/~jinzi/research/convectivePlate-np/Movie.

## References

  * Ahlers et al. (2009) Ahlers, G., Grossmann, S. & Lohse, D. 2009 Heat transfer and large scale dynamics in turbulent Rayleigh-Bénard convection. Rev. Mod. Phys. 81 (2), 503.
  * Araujo et al. (2005) Araujo, F. F., Grossmann, S. & Lohse, D. 2005 Wind Reversals in Turbulent Rayleigh-Bénard Convection. Phys. Rev. Lett. 95 (8), 084502.
  * Bao et al. (2015) Bao, Y., Chen, J., Liu, B.-F., She, Z.-S., Zhang, J. & Zhou, Q. 2015 Enhanced heat transport in partitioned thermal convection. Journal of Fluid Mechanics 784, R5.
  * Brown & Ahlers (2007) Brown, E. & Ahlers, G. 2007 Large-Scale Circulation Model for Turbulent Rayleigh-Bénard Convection. Phys. Rev. Lett. 98 (13), 134501\.
  * Burke (2011) Burke, K. 2011 Plate Tectonics, the Wilson Cycle, and Mantle Plumes: Geodynamics from the Top. Annual Review of Earth and Planetary Sciences 39 (Volume 39, 2011), 1–29.
  * Chernih et al. (2014) Chernih, A., Sloan, I. H. & Womersley, R. S. 2014 Wendland functions with increasing smoothness converge to a Gaussian. Advances in Computational Mathematics 40 (1), 185–200.
  * Elder (1967) Elder, J. 1967 Convective self-propulsion of continents. Nature 214 (5089), 657–660.
  * Gurnis (1988) Gurnis, M. 1988 Large-scale mantle convection and the aggregation and dispersal of supercontinents. Nature 332 (6166), 695–699.
  * Hou & Li (2007) Hou, T. Y. & Li, R. 2007 Computing nearly singular solutions using pseudo-spectral methods. Journal of Computational Physics 226 (1), 379–397.
  * Howard et al. (1970) Howard, L. N., Malkus, W. V. R. & Whitehead, J. A. 1970 Self-convection of floating heat sources: A model for continental drift. Geophysical and Astrophysical Fluid Dynamics 1 (1-2), 123–142.
  * Huang (2024) Huang, J. M. 2024 Covering convection with a thermal blanket: numerical simulation and stochastic modelling. Journal of Fluid Mechanics 980, A47.
  * Huang & Moore (2022) Huang, J. M. & Moore, N. J. 2022 Morphological attractors in natural convective dissolution. Physical Review Letters 128 (2), 024501\.
  * Huang et al. (2021) Huang, J. M., Shelley, M. J. & Stein, D. B. 2021 A stable and accurate scheme for solving the Stefan problem coupled with natural convection using the Immersed Boundary Smooth Extension method. Journal of Computational Physics 432, 110162.
  * Huang & Zhang (2022) Huang, J. M. & Zhang, J. 2022 Rayleigh–Bénard thermal convection perturbed by a horizontal heat flux. J. Fluid Mech. 954.
  * Huang et al. (2018) Huang, J. M., Zhong, J.-Q., Zhang, J. & Mertz, L. 2018 Stochastic dynamics of fluid–structure interaction in turbulent thermal convection. J. Fluid Mech. 854.
  * Kious & Tilling (1996) Kious, W. J. & Tilling, R. I. 1996 This dynamic Earth: the story of plate tectonics. DIANE Publishing.
  * Li et al. (2024) Li, Y.-Z., Chen, X. & Xi, H.-D. 2024 Enhanced heat transfer and reduced flow reversals in turbulent thermal convection with an obstructed centre. Journal of Fluid Mechanics 981, A16.
  * Lowenstein (2024) Lowenstein, J. H. 2024 Near-critical behavior of the Zhong-Zhang model. Phys. Rev. E 109, 025102.
  * Lowman & Gable (1999) Lowman, J. P. & Gable, C. W. 1999 Thermal evolution of the mantle following continental aggregation in 3D convection models. Geophysical Research Letters 26 (17), 2649–2652.
  * Lowman & Jarvis (1993) Lowman, J. P. & Jarvis, G. T. 1993 Mantle convection flow reversals due to continental collisions. Geophysical Research Letters 20 (19), 2087–2090.
  * Lowman & Jarvis (1995) Lowman, J. P. & Jarvis, G. T. 1995 Mantle convection models of continental collision and breakup incorporating finite thickness plates. Physics of the Earth and Planetary Interiors 88 (1), 53–68.
  * Lowman & Jarvis (1999) Lowman, J. P. & Jarvis, G. T. 1999 Effects of mantle heat source distribution on supercontinent stability. Journal of Geophysical Research: Solid Earth 104 (B6), 12733–12746.
  * Lowman et al. (2001) Lowman, J. P., King, S. D. & Gable, C. W. 2001 The influence of tectonic plates on mantle convection patterns, temperature and heat flow. Geophysical Journal International 146 (3), 619–636.
  * Mao (2021) Mao, Y. 2021 An insulating plate drifting over a thermally convecting fluid: the effect of plate size on plate motion, coupling modes and flow structure. Journal of Fluid Mechanics 916, A18.
  * Mao et al. (2019) Mao, Y., Zhong, J.-Q. & Zhang, J. 2019 The dynamics of an insulating plate over a thermally convecting fluid and its implication for continent movement over convective mantle. Journal of Fluid Mechanics 868, 286–315.
  * Mercier et al. (2014) Mercier, M. J., Ardekani, A. M., Allshouse, M. R., Doyle, B. & Peacock, T. 2014 Self-propulsion of immersed objects via natural bioconvection. Physical Review Letters 112 (20), 204501.
  * Meyers et al. (1987) Meyers, R. A. et al. 1987 Encyclopedia of physical science and technology. Academic Press.
  * Moore & Huang (2023) Moore, N. J. & Huang, J. M. 2023 Fluid pendulum explains reversals of the large-scale circulation in thermal convection. arXiv preprint arXiv:2307.13148 .
  * Peyret (2002) Peyret, R. 2002 Spectral methods for incompressible viscous flow, , vol. 148. Springer Science & Business Media.
  * Plummer et al. (2001) Plummer, C. C., McGeary, D., Carlson, D. H. et al. 2001 Physical geology. McGraw-Hill Boston.
  * Rozel et al. (2017) Rozel, A. B., Golabek, G. J., Jain, C., Tackley, P. J. & Gerya, T. 2017 Continental crust formation on early Earth controlled by intrusive magmatism. Nature 545 (7654), 332–335.
  * Trefethen (2000) Trefethen, L. N. 2000 Spectral methods in MATLAB. SIAM.
  * Turcotte & Schubert (2002) Turcotte, D. L. & Schubert, G. 2002 Geodynamics. Cambridge University Press.
  * Wang & Zhang (2023) Wang, K. & Zhang, J. 2023 Persistent corotation of the large-scale flow of thermal convection and an immersed free body. Proceedings of the National Academy of Sciences 120 (21), e2217705120.
  * Weady et al. (2022) Weady, S., Tong, J., Zidovska, A. & Ristroph, L. 2022 Anomalous Convective Flows Carve Pinnacles and Scallops in Melting Ice. Phys. Rev. Lett. 128, 044502.
  * Whitehead (1972) Whitehead, J. A. 1972 Moving heaters as a model of continental drift. Phys. Earth Planet. In. 5, 199–212.
  * Whitehead (2023) Whitehead, J. A. 2023 Convection cells with accumulating crust: models of continent and mantle evolution. Journal of Geophysical Research: Solid Earth 128 (4), e2022JB025643.
  * Whitehead & Behn (2015) Whitehead, J. A. & Behn, M. D. 2015 The continental drift convection cell. Geophys. Res. Lett. 42 (11), 4301–4308.
  * Yang & Song (2023) Yang, Y. & Song, X. 2023 Multidecadal variation of the Earth’s inner-core rotation. Nature Geoscience 16 (2), 182–187.
  * Zhang & Libchaber (2000) Zhang, J. & Libchaber, A. 2000 Periodic boundary motion in thermal turbulence. Physical Review Letters 84 (19), 4361.
  * Zhang et al. (2020) Zhang, S., Xia, Z., Zhou, Q. & Chen, S. 2020 Controlling flow reversal in two-dimensional Rayleigh–Bénard convection. Journal of Fluid Mechanics 891, R4.
  * Zhong & Zhang (2005) Zhong, J.-Q. & Zhang, J. 2005 Thermal convection with a freely moving top boundary. Phys. Fluids 17 (11), 115105.
  * Zhong & Zhang (2007a) Zhong, J.-Q. & Zhang, J. 2007a Dynamical states of a mobile heat blanket on a thermally convecting fluid. Physical Review E 75 (5), 055301.
  * Zhong & Zhang (2007b) Zhong, J.-Q. & Zhang, J. 2007b Modeling the dynamics of a free boundary on turbulent thermal convection. Physical Review E 76 (1), 016307.
  * Zhong & Gurnis (1993) Zhong, S. & Gurnis, M. 1993 Dynamic feedback between a continentlike raft and thermal convection. Journal of Geophysical Research: Solid Earth 98 (B7), 12219–12232.
  * Zhong et al. (2000) Zhong, S., Zuber, M. T., Moresi, L. & Gurnis, M. 2000 Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection. Journal of Geophysical Research: Solid Earth 105 (B5), 11063–11082.


# An Analytics of Culture: Modeling Subjectivity, Scalability, Contextuality,
and Temporality

Nanne van Noord
Informatics Institute
University of Amsterdam
<EMAIL_ADDRESS>
Melvin Wevers11footnotemark: 1
Department of History
University of Amsterdam
<EMAIL_ADDRESS>
Tobias Blanke
Institute for Logic, Language and Computation
University of Amsterdam
<EMAIL_ADDRESS>
Julia Noordegraaf
Department of Media Studies
University of Amsterdam
<EMAIL_ADDRESS>
Marcel Worring
Informatics Institute
University of Amsterdam
<EMAIL_ADDRESS>Joint first authorship. Order of first authors may swap
depending on context.

##  1 Introduction

There is a bidirectional relationship between culture and AI; AI models are
increasingly used to analyse culture, thereby shaping our understanding of
culture. On the other hand, the models are trained on collections of cultural
artifacts thereby implicitly, and not always correctly, encoding expressions
of culture. This creates a tension that both limits the use of AI for
analysing culture and leads to problems in AI with respect to cultural complex
issues such as bias.

One approach to overcome this tension is to more extensively take into account
the intricacies and complexities of culture. We structure our discussion using
four concepts that guide humanistic inquiry into culture: subjectivity,
scalability, contextuality, and temporality. We focus on these concepts
because they have not yet been sufficiently represented in AI research. We
believe that possible implementations of these aspects into AI research leads
to AI that better captures the complexities of culture. In what follows, we
briefly describe these four concepts and their absence in AI research. For
each concept, we define possible research challenges.

##  2 Subjectivity

In analysing culture, we often want to incorporate and compare viewpoints by
different subjects, providing us with a more holistic perspective on culture.
Subjectivity is inherently relational, meaning that viewpoints gain importance
through their relationship with other viewpoints. These subjectivities are not
represented in a monolithic model that is a simple amalgamation of viewpoints.
Consequently, it is difficult to identify and contrast viewpoints in a model.

Frequency-based AI approaches may not always accurately represent the
prevalence of viewpoints in society, leading to an over-representation of
minority and/or undesirable viewpoints. Attempts have been made to de-bias
models by removing undesired viewpoints or to emphasise desirable viewpoints
[2, 10] Rather than simply removing diversity, the research challenge is to
explicitly encode viewpoints into a model, or into multiple models. The result
could be models representing specific personas; from elites to
underrepresented groups. For example, when training an embedding model on
historical textual data, one could encapsulate social class into the model.

##  3 Scalability

Cultural analysis understands cultural phenomena as inherently multi-scale
through time and space, working though the interactions of particularity and
universality. Digital cultural data can be viewed at different levels,
allowing for ‘close’ and ‘distant reading/viewing’. However, this binary
categorization overlooks the interactions between the cultural phenomena on
different levels. We need AI systems to be able to scale down to the micro-
levels of individual interactions and up to the macro-levels of collective
transformation across varying modalities [4, 5].

When modelling cultural data for AI research, we have to take into
consideration these interactions between scales and the commonalities across
them. Otherwise, we run the risk of conflating culture into one big melting
pot, overlooking differences between individuals and groups as well as
continuities across groups. We need theoretical frameworks to explain large-
scale trends through small-scale processes. This also means that we have to
carefully assess which parts of humanistic inquiry can be automated and scaled
using AI systems and which parts require human interpretation or human-machine
interaction.

##  4 Contextuality

Cultural artifacts obtain their meaning and significance from the context in
which they originate, circulate, and are perceived. Traditionally, a
humanities scholar studies the cultural artifact and its context and learns to
intuit what is meaningful and what is not, thereby advancing their
understanding and insights about the artifact. In addition to the artifact
itself and its metadata, therefore, a fuller understanding of its cultural
value requires data that represents these contexts. In a review of context,
Bradley and Dunlop [3] make the distinction between incidental and meaningful
context. Where the incidental context is that which just happens to be
present, whereas meaningful context may influence the process of extracting
information for obtaining insight. The challenge when analysing culture is to
determine which contextual information is incidental and which is meaningful.

Working with archival data requires knowledge and expertise to adequately
determine what is meaningful. Information offered through metadata can only in
part assist in determining what is meaningful. To be able to make the
distinguishing between incidental and meaningful, we need data that goes
beyond the where, when, and by whom an artifact originated, but also includes
the why [6]. Currently, many models ignore the metadata outright and deal with
the artifact data only. The open question is how can computational models be
derived that not only have sufficient ‘base’ knowledge to contextualise but
that are also able to distinguish between the incidental and meaningful?
Additionally, how can we take into account that not all meaningful context
might be available as data? Situated Data Analysis might offer entry points
for this latter question [7].

##  5 Temporality

Time is crucial for cultural analysis in part because the order of events
influences meaning formation [11]. The meaning of concepts relating to sexual
identity, such as queer, or the visual representation of objects, such as
telephones, cannot be accurately modelled without considering temporal
changes. Currently, in analysing cultural data with AI this temporal dimension
is ignored, thereby disregarding processes related to intertextuality,
intericonicity, and the formation of canons [9]. A cultural expression’s
ability to refer to earlier expressions, or its position in a canon determines
how we view them and what role they play.

Datasets are often temporal snapshots, capturing a single moment in time and
the meaning that concepts had at that point [1]. Currently, having multiple
snapshots allows us to contrast and compare different senses, but it does not
allow us to capture their meaning evolved. To be able to capture how meaning
is constituted through cultural processes, AI models need to encode how
meaning forms over time and how meanings are retained and possibly forgotten
over time [8]. Relying on longitudinal data by itself does not offer the
solution, the challenge is in finding the relationality between temporal data
points and the principles driving how cultural data is produced over time.

##  6 Conclusion

We have described four concepts of the humanistic analysis of culture and the
challenges they present for AI. We contend that there is a fruitful trading
zone where humanistic and AI scholarship can meet and shape the development of
AI. As such, this paper provides an invitation for productive engagement and
collaboration between domains to do justice to the complexity of human
culture.

## References

  * Blanke [2018] Tobias Blanke.  Predicting the past.  _Digital Humanities Quarterly_ , 12(2):1–21, Jul 2018.  ISSN 1938-4122.
  * [2] Tolga Bolukbasi, Kai-Wei Chang, James Y Zou, Venkatesh Saligrama, and Adam T Kalai.  Man is to Computer Programmer as Woman is to Homemaker? Debiasing Word Embeddings.  In D. D. Lee, M. Sugiyama, U. V. Luxburg, I. Guyon, and R. Garnett, editors, _Advances in Neural Information Processing Systems 29_ , pages 4349–4357. Curran Associates, Inc.
  * Bradley and Dunlop [2005] Nicholas A Bradley and Mark D Dunlop.  Towards a multidisciplinary model of context to support context-aware computing.  _Journal of Human-Computer Interaction_ , 20(4):403–446, 2005.
  * [4] Simon Carrignon, Tom Brughmans, and Iza Romanowska.  Tableware trade in the Roman East: Exploring cultural and economic transmission with agent-based modelling and approximate Bayesian computation.  15(11):e0240414.
  * [5] Edward D. Lee, Bryan C. Daniels, Christopher R. Myers, David C. Krakauer, and Jessica C. Flack.  Scaling theory of armed-conflict avalanches.  102(4):042312.
  * Noordegraaf [2015] Julia Noordegraaf.  Crowdsourcing television’s past: The state of knowledge in digital archives.  _TMG Journal for Media History_ , 14(22):108–120, Sep 2015.  ISSN 2213-7653.  doi: 10.18146/tmg.139.
  * Rettberg [2020] Jill Walker Rettberg.  Situated data analysis: a new method for analysing encoded power relationships in social media platforms and apps.  _Humanities and Social Sciences Communications_ , 7(11):1–13, Jun 2020.  ISSN 2662-9992.  doi: 10.1057/s41599-020-0495-3.
  * [8] Troy Tassier.  A model of fads, fashions, and group formation.  9(5):51–61.
  * van Noord [2022] Nanne van Noord.  A survey of computational methods for iconic image analysis.  _Digital Scholarship in the Humanities_ , 2022.  ISSN 2055-7671.  doi: 10.1093/llc/fqac003.
  * [10] Melvin Wevers.  Using Word Embeddings to Examine Gender Bias in Dutch Newspapers, 1950-1990.  In _Proceedings of the 1st International Workshop on Computational Approaches to Historical Language Change_ , pages 92–97. Association for Computational Linguistics.
  * [11] Melvin Wevers, Jan Kostkan, and Kristoffer L. Nielbo.  Event Flow - How Events Shaped the Flow of the News, 1950-1995.  In Maud Ehrmann, Folgert Karsdorp, Melvin Wevers, Tara Lee Andrews, Manuel Burghardt, Mike Kestemont, Enrique Manjavacas, Michael Piotrowski, and Joris van Zundert, editors, _Proceedings of the Conference on Computational Humanities Research 2021_ , volume 2989 of _CEUR Workshop Proceedings_ , pages 62–76. CEUR.