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

(2.1)$$\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}$$

(2.2)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$

(2.3)$$\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$:

(2.4)\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

(2.5)\begin{gather} \theta = 0,\quad u_y=v=0 \quad \mbox{for } y=1\ \mbox{and}\ x\notin P, \end{gather}

(2.6)\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

(2.7)\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

(2.8)\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.

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,

(4.1)\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

(4.2)\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

(4.3)\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,

(4.4)\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

(4.5)\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$,

(4.6)\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 )$:

(4.7)$$\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}$$

(4.8)$$\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

(4.9)\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

(4.10)\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,

(5.1)$$\begin{gather} \bar{\theta}(x,t) = \alpha - \alpha\cos{kr[x-x_m(t)]}, \end{gather}$$

(5.2)$$\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$:

(5.3)$$\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}$$

(5.4)$$\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

(5.5)\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

(5.6)$$\begin{gather} \mbox{tr}(J) \to -\lambda <0, \end{gather}$$

(5.7)$$\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

(5.8)\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

(5.9)\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.