1. Introduction
The interior of Earth has fascinated generations of scientists (Plummer, Carlson & Hammersley Reference Plummer, Carlson and Hammersley2015). Among them, Leonardo da Vinci (1452–1519) was one of the pioneers who noticed the incessant geological movements of our planet, as he observed the presence of marine fossils in the mountains. We now know that the continents of Earth do not stay in place and instead undergo tectonic motions, and thermal convection in Earth's mantle is believed to be the driving force of these motions (Kious & Tilling Reference Kious and Tilling1996).
Thermal convection occurs when uneven temperatures of fluid lead to uneven density and buoyancy, so warm fluid rises while cold fluid sinks. The definition of fluids here can be very broad, as modern geologists confirm that even the mantle flows like fluids at a large time scale (Turcotte & Schubert Reference Turcotte and Schubert2002). The Prandtl number (
${\textit {Pr}}{}$) there, defined as the ratio between the mantle's kinematic viscosity and thermal diffusivity, is estimated to be approximately 
$10^{23}$ (Meyers et al. Reference Meyers1987).
The core of Earth is much warmer than its surface, and the destabilizing buoyancy is strong enough to drive mantle convection. As a measure of relative strength between buoyancy and viscous effects, the Rayleigh number (
${\textit {Ra}}{}$) is approximately 
$10^6$ in the mantle (Selley, Cocks & Plimer Reference Selley, Cocks and Plimer2005). In the well-studied case of Rayleigh–Bénard convection, such a high 
${\textit {Ra}}{}$ is known to lead to turbulent fluid motions (Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009). With the mantle convecting like a fluid, its surface flow transports the continental plates, resulting in their tectonic motions.
Due to the large spatial scale of Earth and the long time scale of mantle convection, the geophysical study of plate tectonics focuses on the current state of continents as well as predicting its important consequences, such as earthquakes (Plummer et al. Reference Plummer, Carlson and Hammersley2015). On the other hand, numerical simulations (Howard, Malkus & Whitehead Reference Howard, Malkus and Whitehead1970; Whitehead Reference Whitehead1972, Reference Whitehead2023; Whitehead & Behn Reference Whitehead and Behn2015; Mao, Zhong & Zhang Reference Mao, Zhong and Zhang2019; Mao Reference Mao2021) and lab-scale experiments (Elder Reference Elder1967; Zhang & Libchaber Reference Zhang and Libchaber2000; Zhong & Zhang Reference Zhong and Zhang2005, Reference Zhong and Zhang2007a,Reference Zhong and Zhangb; Whitehead, Shea & Behn Reference Whitehead, Shea and Behn2011) have proven to be effective means of understanding the dynamics of plate tectonics.
Early experiments of Elder (Reference Elder1967) showcase how one can recover the tectonic motion in the lab, where a paraffin fluid layer was used to model the mantle, while a thin sheet of plastic floating on top served as a model continental plate. When the paraffin is heated from below, convection occurs and the plate moves due to shearing of the convective flow beneath. Such a simple experimental set-up also displays non-trivial dynamics; as Zhang & Libchaber (Reference Zhang and Libchaber2000) observed, the plate moves periodically between the two bounding walls of the fluid surface. Through more detailed investigations (Zhong & Zhang Reference Zhong and Zhang2005, Reference Zhong and Zhang2007a,Reference Zhong and Zhangb), the size of the floating plate is shown to strongly affect the plate motion, where small plates display periodic motion, while large plates stay trapped in the middle of the fluid surface. Placing a moving heat source on top of a thermally convecting fluid also yields plate motions, and this is investigated experimentally by Howard et al. (Reference Howard, Malkus and Whitehead1970) and Whitehead (Reference Whitehead1972). Although the geometry, physical parameters and time scales presented in these works are very different from the mantle convection, they reveal surprising dynamics, and most importantly, provide invaluable insight into the fluid–structure interaction mechanism behind continental drift.
The numerical exploration of plate tectonics has developed rapidly in the past several decades. Gurnis (Reference Gurnis1988) provided the first time-dependent numerical simulations of continental drift, where multiple continents were allowed to merge and diverge. As in this work, many other numerical and theoretical endeavours (Zhong & Gurnis Reference Zhong and Gurnis1993; Lowman & Jarvis Reference Lowman and Jarvis1993, Reference Lowman and Jarvis1995, Reference Lowman and Jarvis1999; Lowman & Gable Reference Lowman and Gable1999; Zhong et al. Reference Zhong, Zuber, Moresi and Gurnis2000; Lowman, King & Gable Reference Lowman, King and Gable2001) employ geophysical parameters of the mantle and enable rapid advancement of our understanding of the interior of Earth. It has also become clear that the two-way coupling between the continental plates and the mantle convection results in diverse dynamics (Gurnis Reference Gurnis1988; Zhong & Gurnis Reference Zhong and Gurnis1993; Phillips & Bunge Reference Phillips and Bunge2005; Whitehead & Behn Reference Whitehead and Behn2015). Most notably, large plates are observed to have more consistent motions, while small plates tend to move sporadically (Gurnis Reference Gurnis1988; Whitehead & Behn Reference Whitehead and Behn2015). These observations were examined recently by Mao et al. (Reference Mao, Zhong and Zhang2019) and Mao (Reference Mao2021) through resolved numerical simulations.
The aforementioned works confirm that continental plates are not only passive to the mantle flow advection beneath, but also affecting the flow structure through the thermal blanket effect: The continental crust is known to have a much lower heat flux compared to the oceanic crust due to its large crust depth (Mao Reference Mao2021), so the continental plates serve essentially as a blanket that prevents heat from escaping and warms up the mantle beneath. The warm and light mantle tends to rise, forming an upward convective flow. As this flow moves towards the surface of Earth, it diverges and creates a fluid forcing beneath the continental plate, transporting the plate away. This is the current understanding of continental drift, and the thermal blanket effect has been verified both numerically and experimentally (Gurnis Reference Gurnis1988; Zhong & Zhang Reference Zhong and Zhang2005, Reference Zhong and Zhang2007a,Reference Zhong and Zhangb).
This paper aims to provide a new angle for modelling plate tectonics through a low-dimensional model. After conducting direct numerical simulations (DNS) of the plate-flow interaction in a two-dimensional periodic domain, we propose a stochastic model of the plate motion, and show how the moving plate couples mechanically and thermally to a convecting fluid flow beneath it. Retaining only the most basic physics of thermal convection, this model recovers the dynamics observed in the full DNS, and captures the transition of plate dynamics seen in Gurnis (Reference Gurnis1988), Whitehead & Behn (Reference Whitehead and Behn2015), Mao et al. (Reference Mao, Zhong and Zhang2019) and Mao (Reference Mao2021).
In what follows, we summarize the equations and numerical methods in § 2, and present the numerical results in § 3. The stochastic model is derived systematically in § 4, and its application to the convection domain with various aspect ratios is discussed in § 5. Finally, we summarize and discuss our results in § 6.
2. Numerical model
2.1. Flow equations
The configuration of our numerical simulation is shown in figure 1, where a solid plate centred at location 
$x=x_p$ floats on top of a convecting fluid that is bounded in the 
$y$ direction and periodic in the 
$x$ direction. Throughout this study, all lengths are rescaled by the fluid depth 
$H$, time is rescaled by the diffusion time 
$H^2/\kappa$ (where 
$\kappa$ is thermal diffusivity), and temperature is rescaled by the temperature difference 
$\Delta T$ between the bottom and top free surfaces. The 
$x$ direction of the fluid domain is periodic with period 
$\varGamma = D/H$ (where 
$D$ is the domain width), so the overall computational domain is 
$x\in (0,\varGamma )$ and 
$y\in (0,1)$, as shown in figure 1. With the Boussinesq approximation, the resulting partial differential equations for flow speed 
$\boldsymbol {u} = (u,v)$, pressure 
$p$, and temperature 
$\theta \in [0,1]$ are
$$\begin{gather} \frac{{\rm D}\boldsymbol{u}}{{\rm D}t} ={-}\boldsymbol{\nabla} p + {\textit {Pr}}{}\,{\nabla}^2\boldsymbol{u} + {\textit {Ra}}{}\,{\textit {Pr}}{}\,\theta\,\boldsymbol{e}_y, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$
$$\begin{gather}\frac{{\rm D} \theta}{{\rm D}t} = \nabla^2 \theta. \end{gather}$$
Here, the Rayleigh number is 
${\textit {Ra}}{} = \alpha g\,\Delta T\,H^3 /\nu \kappa$ and the Prandtl number is 
${\textit {Pr}}{} = \nu /\kappa$, where 
$\nu$, 
$\alpha$ and 
$g$ are the kinematic viscosity, the thermal expansion coefficient of the fluid, and the acceleration due to gravity, respectively. Simple modifications to the flow solver can be made for the geophysical mantle convection, but as we wish to consider a more general case of fluid–structure interactions and to apply our theory to future laboratory experiments, we preserve the inertia of both the fluid and the solid plate in this study.
Figure 1. Schematics of the moving plate and convecting fluid. The fluid domain is bounded by 
$y\in (0,1)$ and periodic in 
$x$, and the floating plate of width 
$d$ has centre location 
$x_p$.
2.2. Boundary conditions
Without the presence of a plate, the boundary conditions are straightforward: the flow velocity 
$ {\boldsymbol {u}} = (u,v)$ is no-slip at the fluid/solid boundary, and shear-free at the air/fluid interface. The temperature 
$\theta$ is 1 at the bottom and 0 at the air/fluid boundary.
This yields the boundary conditions for the bottom surface 
$y=0$:
\begin{equation} \theta = 1,\quad u=v=0 \quad \mbox{at } y=0. \end{equation} At the top surface, the plate is effectively shielding the heat from escaping, so we take 
$\theta _n =0$ there, while setting the flow to be no-slip with respect to the moving plate. The resulting boundary conditions are
\begin{gather} \theta = 0,\quad u_y=v=0 \quad \mbox{for } y=1\ \mbox{and}\ x\notin P, \end{gather}
\begin{gather}\theta_y = 0,\quad u = u_p,\quad v=0 \quad \mbox{for } y=1\ \mbox{and}\ x\in P. \end{gather}Alternatively, these conditions can be enforced as
\begin{equation} \left.\begin{array}{ll@{}} (1-\mathbb{1}_{P})\theta + \mathbb{1}_{P}\theta_y = 0,\\ \mathbb{1}_{P}u+(1-\mathbb{1}_{P})u_{y} = u_p & \mbox{at } y=1, \\ v = 0, \end{array}\right\} \end{equation}
where 
$\mathbb {1}_{P}$ is an indicator function that takes the value 1 under the plate and 0 otherwise.
2.3. Plate dynamics
The fluid shear force drives the plate motion directly, so
\begin{equation} m \dot{u}_p ={-}{\textit {Pr}}{} \int_P \frac{ \partial{u}}{ \partial y } (x,1,t)\,{{\rm d}\kern 0.06em x}. \end{equation}
Here, 
$u_p = \dot {x}_p$ is the plate velocity, 
$m = \rho d$ is the dimensionless mass of the plate with linear density 
$\rho$ and width 
$d$, and the integration area 
$P = \{x\, |\, x\in (x_p-d/2, x_p+d/2)\}$ is the region under the plate.
2.4. Parameters and numerical method
The numerical method solving (2.1)–(2.3) with (2.4), (2.7) and (2.8) is detailed in Appendix A, where we use a Fourier–Chebyshev spectral method to obtain resolved and accurate numerical solutions. In all simulations, we choose 
${\textit {Ra}}{} = 10^6$, 
${\textit {Pr}} = 7.9$, 
$\varGamma = 1\unicode{x2013}16$ and 
$m = 4d$ (where 
$d$ is the plate width), matching the parameters of water convection in experiments (Zhang & Libchaber Reference Zhang and Libchaber2000; Zhong & Zhang Reference Zhong and Zhang2005).
3. Numerical results
Several trajectories of plates with various sizes are shown in figure 2(a), and we can see immediately how the plate's size affects its motion. We define the covering ratio 
${\textit {Cr}}{} = d/\varGamma$ to measure how much of the free surface is covered by the plate of width 
$d$, where the fluid aspect ratio is 
$\varGamma = 4$ for all the results in this section. For small plates, their net displacement is small, which can be seen better from the total displacement 
$d_p(t) = \int _0^t |u_p(t')|\,{\rm d}t'$ shown in figure 2(b). Increasing the plate size, linear motion appears as 
${\textit {Cr}}{}$ becomes greater than 0.33, as seen in the green trajectories in figure 2(a). These trajectories are subject to reversals, as there is an effective noise from the turbulent fluid forcing. As 
${\textit {Cr}}{}$ increases further, the linear motion becomes more persistent, as the reversals of plate motion become rare when 
${\textit {Cr}}{}\to 1$ in figure 2(a). We note that similar dynamical behaviours have been seen in geophysical Stokes flow simulations (Mao Reference Mao2021), therefore the coupling mechanism between the moving plate and flow beneath must be similar for different flow regimes.
Figure 2. Motion of a plate floating on top of a convective fluid. (a) Trajectories of plates of various sizes. The covering ratio is 
${\textit {Cr}}{} = d/\varGamma$, and the time 
$t_0$ marks the beginning of plate motion. (b) The total displacement of the plates, where 
$d_p$ is defined through 
$\dot {d}_p = |\dot {x}_p|$. (c) Average plate speed 
$v_p = \langle |\dot {x}_p|\rangle$ becomes high when 
${\textit {Cr}}>0.33$ and reaches a maximum at 
${\textit {Cr}}{} = 0.56$. In these simulations, 
$\varGamma = 4$, 
$\rho =4$, 
${\textit {Ra}}{} = 10^6$, and 
${\textit {Pr}}{} = 7.9$.
From the total displacement 
$d_p$, one can see that a maximum plate speed is achieved at 
${\textit {Cr}}{}\approx 0.5$, and this can be confirmed by plotting the time-averaged plate speed 
$v_p = \langle |\dot {x}_p|\rangle$ in figure 2(c). The average velocity 
$v_p$ remains low for small plates, but increases significantly for 
${\textit {Cr}}{}>0.33$ and reaches a maximum at approximately 
${\textit {Cr}} = 0.5$.
To investigate the transition between dynamical states, the typical flow and temperature distributions in the fluid are shown in figure 3. In figure 3(a), a small plate with 
${\textit {Cr}}{} = 0.125$ is placed on the convecting fluid and is attracted by the centre of downwelling fluid at 
$x_m$ (figure 1), where the surface flow forms a sink. This sink is a stable equilibrium for the plate, as any deviations from this sink will result in a restoring fluid force acting on the plate. The structure of this flow sink can be seen further in figure 3(b), where both the 
$y$-averaged temperature 
$\bar {\theta }= \int _0^1\theta \, {{\rm d}y}$ and the 
$y$-averaged vertical flow velocity 
$\bar {v} = \int _0^1 v\, {{\rm d}y}$ reach their minima.
Figure 3. Dynamics of a floating plate with (a–d) 
${\textit {Cr}}{} = 0.125$ and (e–h) 
${\textit {Cr}}{} = 0.5$. (a) A typical snapshot of the flow and temperature distributions under the plate with 
${\textit {Cr}}{} = 0.125$. (b) The 
$y$-averaged temperature and vertical flow velocity corresponding to (a), the shaded area indicating the position of the plate. (c) The plate displacement is a random process. (d) The plate velocity is a random variable with mean 0, whose probability density function is a Gaussian distribution, as indicated by the histogram. (e) At 
${\textit {Cr}}{} = 0.5$, typical flow and temperature distributions showing that the plate is transported by the surface flow. (f) The 
$y$-averaged temperature and vertical flow velocity corresponding to (e). (g) The plate displacement is linear in time, indicating a unidirectional translation. (h) Plate velocity has non-zero mean. Supplementary movies of these simulations are available at https://doi.org/10.1017/jfm.2023.1071.
Following this surface flow pattern, the plate displacement 
$x_p$ is stochastic, as shown in figure 3(c). Due to the random forcing from turbulent flows, the plate location is subject to noise that can be seen affecting the plate velocity 
$u_p$ in figure 3(d), whose histogram shows a Gaussian distribution. It is rare but not impossible for the plate to experience a strong ‘wind’ from the flow, which can push the plate away from the flow sink, across the flow source, and back to the sink again, resulting in the jumps in 
$x_p$ seen in figure 3(c).
Figure 3(e) shows the dynamics of a plate with 
${\textit {Cr}}{} = 0.5$. In this case, the plate motion is unidirectional, as shown in figures 3(g) and 3(h), with velocity 
$u_p$ that has a non-zero mean. Shown in figures 3(e) and 3(f), the moving plate tends to situate between the flow sink and source. As the surface flow pushes the plate towards its sink, the distribution of flow temperature also shifts, leading to a moving plate chasing a moving surface flow sink.
This is a direct consequence of the thermal blanket effect: When the plate is large enough, the temperature increases beneath it as heat cannot escape there. This local warming modifies the flow temperature and effectively pushes the cold, downwelling fluid away, resulting in a shift of the flow sink location. Overall, the plate moves towards the cold flow sink while simultaneously pushing the sink away. Thus a simple dynamics exists for this seemingly complicated fluid–structure interaction problem, and we will derive a model from these observations.
4. Stochastic model
As seen in figure 3, the variation of the 
$y$-averaged temperature 
$\bar {\theta }$ strongly affects the flow pattern. To capture the variational and periodic nature of 
$\bar {\theta }$, we approximate it with its lowest non-trivial Fourier mode,
\begin{equation} \bar{\theta}(x,t) = \alpha - \alpha\cos{r[x-x_m(t)]}, \end{equation}
where 
$x_m$ is the location of the surface flow sink in figure 1, 
$r = 2{\rm \pi} \varGamma ^{-1}$ is the wavenumber, and the constant 
$\alpha$ measures the strength of temperature variation.
Induced by this temperature distribution, the surface flow velocity 
$U(x,t) = u(x, 1, t)$ can be approximated as
\begin{equation} U(x,t) ={-}\beta\sin{r[x-x_m(t)]}, \end{equation}
where 
$\beta >0$ is the surface flow strength. Indeed, this surface flow profile has a sink at 
$x = x_m$: small deviations from 
$x_m$ result in 
$U>0$ for 
$x< x_m$, and 
$U<0$ for 
$x>x_m$, so locally, the flow points towards 
$x = x_m$.
We note that this surface flow profile does not match the plate velocity at the solid/fluid boundary, and the mismatch between 
$U$ and 
$u_p$ allows us to estimate 
$\partial u/\partial y(x,1,t)$ and the resulting shear stress, leading to the plate acceleration
\begin{equation} \dot{u}_p(t) ={-}\frac{{\textit {Pr}}}{m} \int_P \frac{u_p(t)-U(x,t)}{\delta}\, {{\rm d}\kern 0.06em x} + \sigma\,\dot{W}(t). \end{equation}
Here, 
$\delta$ is the momentum boundary layer thickness (Schlichting & Gersten Reference Schlichting and Gersten2016) that is determined by 
${\textit {Ra}}{}$ and 
${\textit {Pr}}{}$, so we assume it to be constant in our study. We also include a white noise with standard deviation 
$\sigma$, representing the turbulent fluid forcing. Overall, (4.2) and (4.3) represent a flow profile that has both the large-scale circulations and the small-scale turbulent flows, consistent with observations of high-
${\textit {Ra}}{}$ thermal convection (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009; Huang & Zhang Reference Huang and Zhang2022).
To model the moving plate as a thermal blanket, we look at the 
$y$-averaged heat equation,
\begin{equation} \frac{ \partial{\bar{\theta}}}{ \partial t } = \frac{ \partial^2 \bar{\theta}}{ {\partial{x}}^2 } + q(x,t). \end{equation}
Here, we have ignored the flow advection, and 
$q(x,t) = ({\partial \theta }/{\partial y})(x,1,t)-(\partial {\theta }/{\partial y})(x,0,t)$ is the vertical heat flux passing through location 
$x$. Assuming that the heat leaving the fluid–air interface obeys Newton's law of cooling, and no heat penetrates the plate, we can rewrite the heat equation as
\begin{equation} \frac{ \partial{\bar{\theta}}}{ \partial t } = \frac{ \partial^2 \bar{\theta}}{ {\partial{x}}^2 } - \gamma \bar{\theta}(1-\mathbb{1}_P). \end{equation}
The indicator function 
$\mathbb {1}_P(x)$ is 1 when 
$x\in P$ and 0 otherwise, and the constant 
$\gamma$ models the rate of cooling. We now plug in the value of 
$\bar {\theta }$ from (4.1) and integrate over 
$x$, which leads to an ordinary differential equation for 
$x_m$,
\begin{equation} \dot{x}_m(t) = \frac{\gamma}{\rm \pi} \int_P[1-\cos r(x_m-x)]\sin{r(x_m-x)}\, {{\rm d}\kern 0.06em x}. \end{equation} The integrals in (4.3) and (4.6) can be evaluated exactly. Defining a phase angle 
$\phi = r(x_p-x_m)$, we arrive at a closed dynamical system for 
$(u_p, \phi )$:
$$\begin{gather} \dot{u}_p ={-}\frac{\beta\lambda}{{\rm \pi}\,{\textit {Cr}}}\sin({\rm \pi}\,{\textit {Cr}}) \sin\phi - \lambda u_p + \sigma \dot{W}, \end{gather}$$
$$\begin{gather}\dot{\phi} = ru_p + \frac{2\gamma}{\rm \pi}\sin({\rm \pi}\,{\textit {Cr}}) \sin\phi-\frac{\gamma}{2{\rm \pi}}\sin(2{\rm \pi}\,{\textit {Cr}}) \sin(2\phi), \end{gather}$$
where 
$\lambda = {\textit {Pr}}/(\rho \delta )$. Once the dynamics of 
$(u_p, \phi )$ is known, the dynamics of 
$(x_p, x_m)$ can be calculated through 
$\dot {x}_p = u_p$ and 
$x_m = x_p-r^{-1}\phi$.
There are four parameters in this model: 
$\beta$ is the strength of surface flow, 
$\lambda ={\textit {Pr}}/(\rho \delta )$ is a damping coefficient, 
$\gamma$ is the rate of surface cooling, and 
$\sigma$ is the random fluid forcing. Physically, the surface cooling rate 
$\gamma$ is affected by the surface flow strength 
$\beta$, so we take 
$\gamma = r \beta$ in this model, which results in the correct dynamics. The remaining parameters can be estimated from the numerical simulations, and their values and estimation procedures are included in Appendix B.
Starting with 
$\phi (0) = 0$ and a random value of 
$u_p(0)$, figures 4(a)–4(c) show some typical trajectories of 
$x_p(t)$ at different 
${\textit {Cr}}{}$. In figure 4(a), the trajectory of the small plate (
${\textit {Cr}}{} = 0.2$) is noise-driven. For the medium plate with 
${\textit {Cr}}{} = 0.5$, figure 4(b) shows that its trajectory is composed of linear translations with reversals. For the large plate with 
${\textit {Cr}}{} = 0.8$, figure 4(c) indicates that the translation is unidirectional. The plate speed in figure 4(c) is comparable to the full DNS result in figure 2, and this speed decreases as 
${\textit {Cr}}$ increases further. The typical displacement 
$x_p$ for plates with various sizes is shown in figure 4(d), which resembles figure 2(a) and has a transition between the noise-driven and linear motions at 
${\textit {Cr}}{}\approx 0.3$. Thus this simple model captures all the key features of the full numerical simulation.
Figure 4. Simulated trajectories from the stochastic model. (a) The dynamics of the small plate is a random walk. (b) The medium-sized plate has a non-zero translational velocity whose direction is subject to reversals. (c) The translational motion of the large plate is more persistent. (d) Trajectories of all simulated paths, whose dynamics recover figure 2(a).
Without noise, the critical behaviour of the dynamical system (4.7) and (4.8) can be analysed further. For small 
${\textit {Cr}}{}$, it is easy to see that (4.7) and (4.8) have 
$u_p =0$, 
$\phi = 2{\rm \pi} N$ (where 
$N$ is an integer) as equilibria, which are stable and reflect the passive state of plate motion. Increasing 
${\textit {Cr}}{}$, new equilibria appear at 
${\textit {Cr}}{}^* = 1/3$. For 
${\textit {Cr}}{}>{\textit {Cr}}^*$, it becomes possible to have a non-zero plate velocity 
$u_p^* = (\beta /{\rm \pi} )\hat {u}_p^*$, where
\begin{equation} \hat{u}^*_p ={-}\frac{\sin{\rm \pi}\,{\textit {Cr}}}{{\textit {Cr}}}\sin\phi^*, \end{equation}
and the equilibrium 
$\phi ^*$ can be determined from
\begin{equation} \cos{\phi^*} = \frac{1-(2\,{\textit {Cr}})^{{-}1}}{\cos{\rm \pi}\,{\textit {Cr}}}. \end{equation}
These states represent the translation of the plate. We note that (4.9) and (4.10) are functions of only 
${\textit {Cr}}{}$, and are thereby independent of all other parameters assumed in this model. To recover the dimensional plate velocity 
$u_p^*$, one only needs to know the flow speed factor 
$\beta /{\rm \pi}$.
The possible values of 
$\hat {u}_p^*$ and 
$\phi ^*$ are shown in figure 5. For 
${\textit {Cr}}{}<{\textit {Cr}}^* = 1/3$, 
$u^*_p = 0$ and 
$\phi ^* = 2{\rm \pi} N$ are the only possible equilibria that reflect the passive nature of the small plate that is always attracted by the surface flow sink. For 
${\textit {Cr}}{}>{\textit {Cr}}^*$, new phases appear as 
$\phi ^* = (2N+1){\rm \pi} \pm \arccos {([(2\,{\textit {Cr}})^{-1}-1](\cos {\rm \pi}\,{\textit {Cr}})^{-1})}$, which are solutions of (4.10) and become stable for large 
${\textit {Cr}}{}$. They indicate that the larger plate tends to sit between the surface flow sink (
$\phi =2N{\rm \pi}$) and source [
$\phi =(2N+1){\rm \pi}$], confirming our observations in figure 3. As the surface flow points from its source to its sink (see arrows in figure 5b), these new phases indicate two possible plate velocities that are given by (4.9) and shown in figure 5(a). Furthermore, figure 5(a) resembles figure 2(c) as the plate velocity vanishes for 
${\textit {Cr}}\to {\textit {Cr}}^*$ and 
${\textit {Cr}}\to 1$, and obtains its maximum at approximately 
${\textit {Cr}} = 0.5$.
Figure 5. Equilibrium values of the plate velocity and phase angle: (a) dimensionless plate velocity 
$\hat {u}_p$ is 0 for small 
${\textit {Cr}}$, but becomes translational as 
${\textit {Cr}}$ increases; (b) the phase angle 
$\phi$. Blue lines represent the passive state where the plate is attracted by the flow sink, and red/orange curves show the equilibrium phase of the translational state. The surface flow direction is labelled with arrows.
Through this simple model, we see clear physics of how the solid plate interacts with the fluid beneath, and the covering ratio 
${\textit {Cr}}{}$ serves as a measure of the strength of the thermal blanket effect. For small 
${\textit {Cr}}{}$, the thermal blanket effect is weak, thus the plate motion is passive to the fluid. Increasing 
${\textit {Cr}}{}$ beyond 
${\textit {Cr}}^*$, the thermal blanket effect is strong enough to alter the flow and temperature distributions, generating an upwelling that can lead to plate motion. Both figures 2(c) and 5(a) suggest that the plate speed peaks at 
${\textit {Cr}}{}\approx 0.5$ and vanishes as 
${\textit {Cr}}{}\to 1$, which also reflects the competition between the thermal blanket effect and the flow convection. As 
${\textit {Cr}}{}$ increases, more area of the free surface is covered by the plate, and the fluid force beneath is averaged in a larger domain. As this domain may cover both the upwelling and downwelling flows, the total fluid force is affected by 
${\textit {Cr}}$. This can be seen in (4.3), where the surface flow velocity 
$U$ provides plate acceleration. In (4.2), 
$U$ is modelled as sinusoidal, and this profile is integrated in (4.3), therefore increasing the integration area in (4.3) to half of the open surface (
${\textit {Cr}} = 1/2$) will cover the highest contribution of 
$U$. Further increasing the covering area will thereby decrease the contribution of 
$U$, as the integration domain is more than half a period of the sine function. In an extreme case, 
${\textit {Cr}} = 1$ indicates an integration of 
$U$ for a full period, leading to zero fluid force, as seen in figure 5(a).
5. Aspect ratio effect
All the results discussed so far focus on the convective fluid domain of aspect ratio 
$\varGamma = 4$, which matches the geometry of many experimental investigations. Varying the domain aspect ratio will certainly affect the dynamics of the convecting fluid and moving plate, as a more complicated, multi-roll flow structure emerges (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009). In this section, we investigate the plate dynamics at 
$\varGamma =2$, 4, 8 and 16, while keeping other dynamical parameters the same as described in §§ 2–4.
Figure 6 shows some typical temperature distributions of DNS results at 
$\varGamma = 2\unicode{x2013}16$. Clearly, a more complicated flow structure emerges as more convection cells appear with increasing 
$\varGamma$. Figure 7 shows trajectories of the plate at different 
$\varGamma$; a common feature re-emerges as the small plate moves little while the large plate translates.
Figure 6. Snapshots of the convective fluid at various aspect ratios: (a) 
$\varGamma = 2$, (b) 
$\varGamma = 8$, (c) 
$\varGamma = 16$. In these simulations, all other parameters are fixed at 
${\textit {Cr}}{} = 0.2$, 
$\rho =4$, 
${\textit {Ra}}{} = 10^6$ and 
${\textit {Pr}}{} = 7.9$.
Figure 7. Simulated trajectories from the DNS at aspect ratios (a) 
$\varGamma = 2$, (b) 
$\varGamma = 8$, and (c) 
$\varGamma = 16$. In these simulations, 
$\rho =4$, 
${\textit {Ra}}{} = 10^6$ and 
${\textit {Pr}}{} = 7.9$, as in figures 2 and 3. (d) The critical 
${\textit {Cr}}^*$ separating the passive and translating states has a 
$-2/3$ power-law scaling with 
$\varGamma$.
To extend our simple model to cases of high aspect ratio, we consider a bulk temperature and surface flow profile with a more complicated spatial dependence,
$$\begin{gather} \bar{\theta}(x,t) = \alpha - \alpha\cos{kr[x-x_m(t)]}, \end{gather}$$
$$\begin{gather}U(x,t) ={-}\beta\sin{kr[x-x_m(t)]}, \end{gather}$$
where the integer 
$k$ is the most dominant wavenumber in the Fourier spectrum. The inclusion of wavenumber 
$k$ is inspired by the fact that the convection might have a multi-roll structure, so the temperature and flow profiles above indicate that there are 
$2k$ convection rolls in the fluid, with 
$k$ surface sinks and 
$k$ surface sources. We track one of the surface sinks 
$x_m(t)$, which is also the location of lowest fluid temperature.
Through similar derivations as in § 4, we can again obtain a stochastic dynamical system for the phase 
$\phi = r(x_p - x_m)$ and the plate velocity 
$u_p$:
$$\begin{gather} \dot{u}_p ={-}\frac{\beta\lambda}{k{\rm \pi}\,{\textit {Cr}}}\sin(k{\rm \pi}\,{\textit {Cr}}) \sin(k\phi) - \lambda u_p + \sigma \dot{W}, \end{gather}$$
$$\begin{gather}\dot{\phi} = ru_p + \frac{2\gamma}{k^2{\rm \pi}}\sin(k{\rm \pi}\,{\textit {Cr}}) \sin(k \phi)-\frac{\gamma}{2k^2{\rm \pi}}\sin(2k{\rm \pi}\,{\textit {Cr}}) \sin(2k \phi). \end{gather}$$ Without noise, it is easy to verify that (5.3) and (5.4) have passive equilibria 
$u_p^* = 0$ and 
$\phi ^* = 2m{\rm \pi} /k$, where 
$m$ is any integer. The Jacobian of such passive states is
\begin{equation} J = \begin{bmatrix} -\lambda & -\beta\lambda({\rm \pi}\,{\textit {Cr}})^{{-}1}\sin(k{\rm \pi}\,{\textit {Cr}}) \\ r & 2\gamma(k{\rm \pi})^{{-}1} \sin(k{\rm \pi}\,{\textit {Cr}})(1-\cos(k{\rm \pi}\,{\textit {Cr}})) \end{bmatrix}. \end{equation} Evaluating the trace and determinant of 
$J$ in the limit of 
${\textit {Cr}}\to 0$ yields
$$\begin{gather} \mbox{tr}(J) \to -\lambda <0, \end{gather}$$
$$\begin{gather}\mbox{det}(J) \to kr\beta\lambda > 0. \end{gather}$$
Therefore both eigenvalues of 
$J$ are negative, indicating stable passive states for small 
${\textit {Cr}}{}$, and confirming our numerical observations.
The stability of passive states is lost when one or both eigenvalues of 
$J$ become positive, therefore 
$\mbox {det}(J) = 0$ provides an equation to determine the critical 
${\textit {Cr}}^*$. After simplification, we have
\begin{equation} \frac{{\textit {Cr}}^*}{k}\,(1-\cos k{\rm \pi}\,{\textit {Cr}}^*) = \frac{\beta r}{2 \gamma}. \end{equation} The root of (5.8) can be determined numerically, given that the parameters 
$\beta$ and 
$\gamma$ are known. In laboratory experiments and geophysical plate tectonics, the surface flow speed scale 
$\beta$ and ventilation coefficient 
$\gamma$ can be estimated properly and used to determine 
${\textit {Cr}}^*$, therefore (5.8) offers the possibility of determining plate mobility through physical parameters.
We note that 
${\textit {Cr}}^*$ is usually small for large 
$\varGamma$, as shown in figure 7, and the root of (5.8) in this limit can be shown as
\begin{equation} {\textit {Cr}}^*\approx\left(\frac{\beta r}{{\rm \pi}^2 \gamma k}\right)^{1/3} \sim \varGamma^{{-}2/3}. \end{equation}
Here, we have used the relations 
$r \sim \varGamma ^{-1}$ and 
$k\sim \varGamma$, with the latter indicating that the number of convection rolls is proportional to the aspect ratio 
$\varGamma$. Equation (5.9) can indeed be verified – as shown in figure 7(d), the critical 
${\textit {Cr}}^*$ measured from DNS data does follow the 
$-2/3$ power law with 
$\varGamma$.
Much more can be investigated through the dynamics of (5.3) and (5.4), and future experimental studies can certainly be used to address further the interaction between the plate and the convective fluid below.
6. Discussion
In this work, we have explored numerically and theoretically the mechanical and thermal coupling between a moving plate and a convecting fluid beneath it. Inspired by present and past works (Gurnis Reference Gurnis1988; Whitehead & Behn Reference Whitehead and Behn2015; Mao et al. Reference Mao, Zhong and Zhang2019; Mao Reference Mao2021), we propose a stochastic model showing that the plate size (covering ratio) is a deciding factor for the strength of the thermal blanket effect. For small plates, the shielding effect is weak and the plate motion is passive to the flow structure; for large plates, the flow beneath becomes warm enough and an upwelling centre is formed, pushing the plate away and resulting in its translation. The proposed simple model consists of minimal assumptions about the flow and its mechanical and thermal coupling to the plate; however, it is capable of predicting the dynamical transition of the plate motion. Although the direct numerical simulations (DNS) are conducted at a parameter regime similar to laboratory experiments (Zhong & Zhang Reference Zhong and Zhang2005, Reference Zhong and Zhang2007a,Reference Zhong and Zhangb), this stochastic model is Reynolds-number-independent, therefore allowing it to be applied to both laboratory and geophysical scenarios.
Laboratory-scale experiments are usually conducted in a bounded convection system, therefore the plate is limited to move between two walls. The flow structure and its coupling to the plate is very different in bounded convection, as previous works show that the flow has two counter-rotating large-scale circulations, whose strengths are modulated by the location of the plate (Zhong & Zhang Reference Zhong and Zhang2005, Reference Zhong and Zhang2007a,Reference Zhong and Zhangb). The plate is also bounded by the solid boundaries, so it stops moving once it touches the wall. Modelling the plate dynamics in this case is not trivial, and advanced tools such as stochastic variational inequalities (Huang et al. Reference Huang, Zhong, Zhang and Mertz2018) can serve as a mathematical means of analysing such plate–boundary interactions. A periodic domain of convecting fluid is necessary to verify experimentally our stochastic model, and such an experiment would also resemble the mantle convection more closely. Future experiments consisting of an annular convection domain could provide more details to validate the stochastic model.
Although here we investigate only the dynamics of a single plate, our numerical method is capable of handling multiple plates provided that their interactions can be modelled properly. We are investigating such interactions currently, which has led to even more diverse and unpredictable dynamics. For example, if two small plates each have 
${\textit {Cr}}{}<1/3$ but their combined size reaches 
${\textit {Cr}}{}>1/3$, then we have seen that each individual plate moves randomly, but their combined ‘super continent’ can translate. In the geophysical case of plate tectonics, plate interactions are the reason for many volcanic activities and mountain formations, therefore understanding the converging and diverging motion of nearby plates might offer new insights into the fluid mechanics behind these geophysical events.
For simplicity, the simulation and model in this work are both two-dimensional, and extending our results to three dimensions is a current priority. Using the Chebyshev–Fourier–Fourier method, we have implemented the numerical solver for evolving a plate sitting on top of a three-dimensional domain that is periodic in two horizontal directions. Moreover, the simulation of plate tectonics on a spherical shell is also possible through the Chebyshev–Chebyshev–Fourier method, which is a configuration closer to the geophysical plate tectonics. Through analysing the DNS results there, we wish to develop our stochastic model further, and use it to address the fluid–structure interactions happening inside Earth.
Finally, we note that the geophysical plate tectonics is much more complicated than any experiments or numerical simulations conducted so far, as the interior of Earth is such a complex environment and is still being explored by modern science. But although current simple models cannot capture fully the dynamics of continental drifts, we hope that they can still offer some fluid mechanical insights into the geophysics of Earth.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2023.1071.
Funding
J.M.H. thanks J.-Q. Zhong, J. Chu and J. Zhang for useful discussions. J.M.H. acknowledges support from the National Natural Science Foundation of China (12272237, 92252204) and the Science and Technology Commission of Shanghai Municipality (21YF1432100).
Declaration of interests
The author reports no conflict of interest.
Appendix A. Numerical method
The Navier–Stokes Boussinesq equations (2.1)–(2.3) in two dimensions can be written in the vorticity and stream function format as
$$\begin{gather} \frac{{\rm D} \omega}{{\rm D}t} = {\textit {Pr}}{}\, \nabla^2\omega + {\textit {Pr}}{}\, {\textit {Ra}}{}\, \frac{\partial \theta}{\partial x}, \end{gather}$$
$$\begin{gather}-\nabla^2 \psi = \omega,\quad \boldsymbol{u} = \boldsymbol{\nabla}_\perp \psi, \end{gather}$$
$$\begin{gather}\frac{{\rm D} \theta}{{\rm D}t} = \nabla^2 \theta, \end{gather}$$
where the 
$z$-component of vorticity 
$\omega = \hat {\boldsymbol {z}}\boldsymbol {\cdot }\boldsymbol {\nabla }\times \boldsymbol {u}$ and the stream function defined by 
$\boldsymbol {u} = \boldsymbol {\nabla }_\perp \psi = (\psi _y, -\psi _x)$ are solved for, alleviating the difficulty of solving for pressure 
$p$.
In the vorticity and stream function format, the boundary conditions (2.4) and (2.7) can be enforced as
\begin{equation} \theta = 1,\quad \psi=\psi_y=0 \quad \mbox{at } y=0 \end{equation}and
\begin{equation} \left.\begin{array}{ll@{}} (1-\mathbb{1}_{P})\,\theta + \mathbb{1}_{P}\,\theta_y = 0,\\ \mathbb{1}_{P}\,\psi_y+(1-\mathbb{1}_{P})\,\psi_{yy} = u_p & \mbox{at } y=1,\\ \psi = 0. \end{array}\right\} \end{equation}
In the numerical simulations, we soften the edge of indicator function 
$\mathbb {1}_{P}$, making it smoother in order to reduce numerical error.
The time derivatives in (A1) and (A3) are approximated with the second-order Adams–Bashforth backward differentiation method. At time step 
$t_n = n\,\Delta T$, we denote 
$\omega _n(x,y) = \omega (x,y,n\,\Delta T)$, 
$\psi _n(x,y) = \psi (x,y,n\,\Delta T)$ and 
$\theta _n(x,y) = \theta (x,y,n\,\Delta T)$, and (A1)–(A3) become
$$\begin{gather} \nabla^2\omega_n -\sigma_1\omega_n = f_n, \end{gather}$$
$$\begin{gather}\nabla^2 \theta_n -\sigma_2 \theta_n = h_n, \end{gather}$$
$$\begin{gather}-\nabla^2 \psi_n = \omega_n, \end{gather}$$where
$$\begin{gather} \sigma_1 = \frac{3}{2\,{\textit {Pr}}\,\Delta t},\quad \sigma_2 = \frac{3}{2\,\Delta t}, \end{gather}$$
$$\begin{gather}f_n = {\textit {Pr}}^{{-}1}\left[ 2 ({\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \omega)_{n-1} - ({\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \omega)_{n-2}\right] - (2\, {\textit {Pr}}\, \Delta t)^{{-}1} \left(4\omega_{n-1}-\omega_{n-2}\right) -{\textit {Ra}} \left( \frac{ \partial{\theta}}{ \partial x } \right)_{n}, \end{gather}$$
$$\begin{gather}h_n = \left[ 2 ({\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \theta)_{n-1} - ({\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \theta)_{n-2}\right] - (2\,\Delta t)^{{-}1} \left(4\theta_{n-1}-\theta_{n-2}\right)\!. \end{gather}$$Equations (A6)–(A8), together with the inhomogeneous Robin boundary conditions (A4) and (A5), are Helmholtz equations that can be solved by standard spectral methods (Peyret Reference Peyret2002). More details of this numerical solver will be included in future publications.
Nonlinear terms such as 
$ {\boldsymbol {u}} \boldsymbol {\cdot } \boldsymbol {\nabla } \theta$ and 
$ {\boldsymbol {u}} \boldsymbol {\cdot } \boldsymbol {\nabla } \omega$ in (A10) and (A11) are computed pseudo-spectrally with a simple and efficient anti-aliasing filter (Hou & Li Reference Hou and Li2007). With given initial and boundary data, (A7) can be solved first to obtain 
$\theta _{n}$, which is inserted in 
$f_{n}$ so (A6) can be solved next. Finally, (A8) is solved with the known 
$\omega _{n}$.
After solving for the flow and temperature fields, the plate acceleration can be determined as
\begin{equation} a_{p,n} ={-}\frac{{\textit {Pr}}{}}{m} \int_P \frac{ \partial^2 \psi_{n}}{ {\partial{y}}^2 } (x,1) \,{{\rm d}\kern 0.06em x}. \end{equation} The plate velocity 
$u_{p,n}$ and plate location 
$x_{p,n}$ can then be computed through a second-order Adams–Bashforth method:
$$\begin{gather} x_{p,n} = x_{p,n-1}+\frac{\Delta t}{2}\,( 3u_{p,n-1} - u_{p,n-2}), \end{gather}$$
$$\begin{gather}u_{p,n} = u_{p,n-1}+\frac{\Delta t}{2}\,( 3a_{p,n-1} - a_{p,n-2}). \end{gather}$$At 
$\varGamma = 4$, there are typically 
$256$ Fourier modes in the 
$x$ direction, and 
$64$ Chebyshev nodes in the 
$y$ direction, with time step size 
$\Delta t = 10^{-6}$. These parameters are tested to yield resolved and accurate numerical solutions.
Appendix B. Model parameters
The four parameters involved in the stochastic model can be estimated from the DNS data and auxiliary numerical tests. The detailed procedures are listed below.
(i) Parameter
$\beta \approx 400$ is measured directly from the numerical simulation in figure 3.(ii) Parameter
$\lambda \approx 200$ is estimated from ${\textit {Pr}} = 7.9$, $\rho = 4$ and $\delta = 0.01$. The boundary layer thickness $\delta$ is estimated from the relation $\delta \sim (2\,\textit {Nu})^{-1}$, where the Nusselt number $\textit {Nu}{}$ is of the order of $10^1$ as measured from the simulation.(iii) Parameter
$\sigma \approx 200$ is estimated from the variance of the plate centre $x_p$ for small ${\textit {Cr}}{}$. From (4.7), we have $\sigma ^2\approx 2\lambda \,\mbox {Var}(x_p)$, where $\mbox {Var}(x_p)\approx 100$ is measured from the numerical simulation shown in figure 3.
We note that 
$\beta$ and 
$\sigma$ model the surface flow, therefore can be estimated from numerical simulations without floating plate or calculated through scaling relations (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009).
Open access
