Covering convection with a thermal blanket: numerical simulation and stochastic modelling

Abstract Adding moving boundaries to convective fluids is known to result in non-trivial and surprising dynamics, leading to spectacular geoformations ranging from kilometre-scale karst terrains to planetary-scale plate tectonics. On the one hand, the moving solid alters the surrounding flow field, but on the other hand, the flow modifies the motion and shape of the solid. This leads to a two-way coupling that is significant in the study of fluid–structure interactions and in the understanding of geomorphologies. In this work, we investigate the coupling between a floating plate and the convective fluid below it. Through numerical experiments, we show that the motion of this plate is driven by the flow beneath. However, the flow structure is also modified by the presence of the plate, leading to the ‘thermal blanket’ effect where the trapped heat beneath the plate results in buoyant and upwelling flows that in turn push the plate away. By analysing this two-way coupling between moving boundary and fluid, we are able to capture the dynamical behaviours of this plate through a low-dimensional stochastic model. Geophysically, the thermal blanket effect is believed to drive the continental drift, therefore understanding this mechanism has significance beyond fluid dynamics.


Introduction
The interior of Earth has fascinated generations of scientists (Plummer et al. 2001).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 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 1996).
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 2002).The Prandtl number (Pr) there, defined as the ratio between the mantle's kinematic viscosity and thermal diffusivity, is estimated to be around 10 23 (Meyers et al. 1987).
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 (Ra) is around 10 6 in the mantle (Selley et al. 2005).In the well-studied case of Rayleigh-Bénard convection, such a high Ra is known to lead to turbulent fluid motions (Ahlers et al. 2009).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 like earthquakes (Plummer et al. 2001).On the other hand, numerical simulations (Howard et al. 1970;Whitehead 1972;Whitehead & Behn 2015;Mao et al. 2019;Mao 2021;Whitehead 2023) and lab-scale experiments (Elder 1967;Zhang & Libchaber 2000;Zhong & Zhang 2005, 2007a,b;Whitehead et al. 2011) have proven to be an effective means of understanding the dynamics of plate tectonics.
Early experiments of Elder (1967) 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 setup also displays nontrivial dynamics, as Zhang & Libchaber (2000) observed the plate moves periodically between the two bounding walls of the fluid surface.Through more detailed investigations (Zhong & Zhang 2005, 2007a,b), the size of 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 experimentally investigated by Howard et al. (1970); Whitehead (1972).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 (1988) 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 endeavors (Zhong & Gurnis 1993;Lowman & Jarvis 1993, 1995;Lowman & Gable 1999;Lowman & Jarvis 1999;Lowman et al. 2001;Zhong et al. 2000) 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 1988;Zhong & Gurnis 1993;Phillips & Bunge 2005;Whitehead & Behn 2015).Most notably, large plates are observed to have more consistent motions, while small plates tend to move sporadically (Gurnis 1988;Whitehead & Behn 2015).These observations were recently examined by Mao et al. (2019); Mao (2021) through resolved numerical simulations.
The aforementioned works confirm that continental plates are not only passive to the mantle flow advection beneath, but are 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 2021), so the continental plates essentially serve 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 1988;Zhong & Zhang 2005, 2007a,b).
This manuscript aims to provide a new angle for modeling plate tectonics through a low-dimensional model.After conducting direct numerical simulations (DNS) of the plate-flow interaction in a 2-dimensional periodic domain, we propose a stochastic model of the plate motion and show how the moving plate mechanically and thermally couples to a convecting fluid flow beneath it.Only retaining 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 (1988); Whitehead & Behn (2015); Mao et al. (2019); Mao (2021).
In what follows, we summarize the equations and numerical methods in Section 2, and present the numerical results in Section 3. The stochastic model will be systematically derived in Section 4, and its application to the convection domain with various aspect ratios will be discussed in Section 5. Finally, we summarize and discuss our results in Section 6.

Flow equations
The configuration of our numerical simulation is shown in Fig. 1, where a solid plate centered 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 /κ (κ is thermal diffusivity), and temperature is rescaled by the temperature difference ∆T between the bottom and top free surface.The x direction of the fluid domain is periodic with period Γ = D/H (D is the domain width), so the overall computational domain is x ∈ (0, Γ ) and y ∈ (0, 1) as shown in Fig. 1.With the Boussinesq approximation, the resulting PDEs for flow speed u = (u, v), pressure p, and temperature θ ∈ [0, 1] are

Du Dt
= −∇p + Pr ∇ 2 u + Ra Pr θ e y , (2.1) Here, the Rayleigh number is Ra = αg∆T H 3 /νκ and the Prandtl number is Pr = ν/κ, where ν, α and g are the kinematic viscosity, the thermal expansion coefficient of the fluid, and the acceleration due to gravity.Simple modifications to the flow solver can be adapted 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.

Boundary conditions
Without the presence of a plate, the boundary conditions are straightforward: the flow velocity u = (u, v) is no-slip at the fluid/solid boundary and shear-free at the air/fluid interface; The temperature θ is 1 at the bottom and 0 at the air/fluid boundary.
This yields the boundary conditions for the bottom surface y = 0, At the top surface, the plate is effectively shielding the heat from escaping, we thus take θ n = 0 there while set the flow to be no-slip with respect to the moving plate.The resulting boundary conditions are θ = 0, u y = v = 0 for y = 1 and x / ∈ P, (2.5) Alternatively, these conditions can be enforced as where 1 P is an indicator function that takes the value of 1 under the plate and 0 otherwise.

Plate dynamics
The fluid shear force directly drives the plate motion, so (2.8) Here u p = ẋp is the plate velocity, m = ρd is the dimensionless mass of the plate with linear density ρ and width d, and the integration area is the region under the plate.

Numerical results
Several trajectories of plates with various sizes are shown in Fig. 2(a), and we can immediately see how the plate's size affects its motion.We define the covering ratio Cr = d/Γ to measure how much of the free surface is covered by the plate of width d, where the fluid aspect ratio is Γ = 4 for all the results in this section.For small plates, their net displacement is small, which can be better seen from the total displacement Increasing the plate size, linear motion appears as Cr becomes greater than 0.33, as seen in the green trajectories in Fig. 2(a).These trajectories are subject to reversals, as there is an effective noise from the turbulent fluid forcing.As Cr further increases, the linear motion becomes more persistent, as the reversals of plate motion become rare when Cr → 1 in Fig. 2(a).We note that similar dynamical behaviors have been seen in geophysical Stokes flow simulations (Mao 2021), therefore the coupling mechanism between the moving plate and flow beneath must be similar for different flow regimes.
From the total displacement d p , one can see that a maximum plate speed is achieved at around Cr ≈ 0.5, and this can be confirmed by plotting the time-averaged plate speed v p = ⟨| ẋp |⟩ in Fig. 2(c).The average velocity v p remains low for small plates, but increases significantly for Cr > 0.33 and reaches a maximum around Cr = 0.5.
To investigate the transition between dynamical states, the typical flow and temperature distributions in the fluid are shown in Fig. 3.In Fig. 3(a), a small plate with Cr = 0.125 is placed on the convecting fluid and it is attracted by the center of downwelling fluid at x m [Fig.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 further seen in Fig. 3(b), where both the y-averaged temperature θ = 1 0 θ dy and the y-averaged vertical flow velocity v = 1 0 v dy reach their minima.Following this surface flow pattern, the plate displacement x p is stochastic as shown in Fig. 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 Fig. 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 Fig. 3(c).
Figure 3(e) shows the dynamics of a plate with Cr = 0.5.In this case, the plate motion is unidirectional as shown in Fig. 3(g)-(h), with velocity u p that has a nonzero mean.Shown in Fig. 3(e)-(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.

Stochastic model
As seen in Fig. 3, the variation of the y-averaged temperature θ strongly affects the flow pattern.To capture the variational and periodic nature of θ, we approximate it with its lowest nontrivial Fourier mode, where x m is the location of surface flow sink in Fig. 1, r = 2πΓ −1 is the wavenumber, and the constant α 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 where β > 0 is the surface flow strength.Indeed, this surface flow profile has a sink at x = x m : Small deviations from x m results in U > 0 for x < x m and U < 0 for x > x m , so the flow locally 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 ∂u/∂y(x, 1, t) and the resulting shear stress, leading to the plate acceleration Here, δ is the momentum boundary layer thickness (Schlichting & Gersten 2016) that is determined by Ra and Pr, so we assume it to be constant in our study.We also include a white noise with standard deviation σ, representing the turbulent fluid forcing.Overall, Eqs.(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-Ra thermal convection (Ahlers et al. 2009;Huang & Zhang 2022).
To model the moving plate as a thermal blanket, we look at the y-averaged heat equation, Here we have ignored the flow advection, and q(x, t) = ∂θ ∂y (x, 1, t) − ∂θ ∂y (x, 0, t) is the vertical heat flux passing through location x.Assuming 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) The indicator function 1 P (x) is 1 when x ∈ P and 0 otherwise, and the constant γ models the rate of cooling.We now plug in the value of θ from Eq. (4.1) and integrate over x, which leads to an ODE for where λ = Pr/(ρδ).Once the dynamics of (u p , ϕ) is known, the dynamics of (x p , x m ) can be calculated through ẋp = u p and x m = x p − r −1 ϕ.
There are four parameters in this model: β as the strength of surface flow, λ = Pr/(ρδ) as a damping coefficient, γ as the rate of surface cooling, and σ as the random fluid forcing.Physically, the surface cooling rate γ is affected by the surface flow strength β, so we take γ = rβ 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 ϕ(0) = 0 and a random value of u p (0), Fig. 4(a)-(c) show some typical trajectories of x p (t) at different Cr .In Fig. 4(a), the trajectory of the small plate (Cr = 0.2) is noise-driven.For the medium plate of Cr = 0.5, Fig. 4(b) shows its trajectory is composed of linear translations with reversals.For the large plate of Cr = 0.8, Fig. 4(c) indicates that the translation is unidirectional.The plate speed in Fig. 4(c) is comparable to the full DNS result in Fig. 2, and this speed decreases as Cr further increases.The typical displacement x p for plates with various sizes is shown in Fig. 4(d), which resembles  Fig. 2(a) and has a transition between the noise-driven and linear motions at Cr ≈ 0.3.Thus, this simple model captures all the key features of the full numerical simulation.Without noise, the critical behavior of the dynamical system Eqs.(4.7) and (4.8) can be further analyzed.For small Cr , it is easy to see that Eqs.(4.7) and (4.8) have u p = 0, ϕ = 2πN (where N is an integer) as equilibria, which are stable and reflect the passive state of plate motion.Increasing Cr , new equilibria appear at Cr * = 1/3.For Cr > Cr  Through this simple model, we see clear physics of how the solid plate interacts with the fluid beneath, and the covering ratio Cr serves as a measure of the strength of the thermal blanket effect.For small Cr , the thermal blanket effect is weak thus the plate motion is passive to the fluid.Increasing Cr beyond 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 Fig. 2(c) and Fig. 5(a) suggest that the plate speed peaks at Cr ≈ 0.5 and vanishes as Cr → 1, which also reflects the competition between the thermal blanket effect and the flow convection.As 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 Cr .This can be seen in Eq. (4.3), where the surface flow velocity U provides plate acceleration.In Eq. (4.2), U is modelled as sinusoidal and this profile is integrated in Eq. (4.3), therefore increasing the integration area in Eq. ( 4.3) to half of the open surface (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, Cr = 1 indicates an integration of U for a full period, leading to zero fluid force as seen in Fig. 5(a).

Aspect ratio effect
All the results discussed so far focus on the convective fluid domain of aspect ratio Γ = 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. 2009).In this section, we investigate the plate dynamics at Γ = 2, 4, 8, and 16 while keeping other dynamical parameters the same as described in Sections 2 to 4.
Figure 6 shows some typical temperature distributions of DNS results at Γ = 2 − 16.Clearly, a more complicated flow structure emerges as more convection cells appear with increasing Γ .Figure 7 shows trajectories of the plate at different Γ ; a common feature reemerges as the small plate moves little while the large plate translates.
To extend our simple model to cases of high aspect ratio, we consider a bulk temper- ature and surface flow profile with a more complicated spatial dependence, (5.1) where the integer k is the most dominant wave number in the Fourier spectrum.The inclusion of wave number 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 Section 4, we can again obtain a stochastic dynamical system for the phase ϕ = r(x p − x m ) and the plate velocity . (5.5) Evaluating the trace and determinant of J in the limit of Cr → 0 yields tr(J) → −λ < 0, (5.6) det(J) → krβλ > 0. (5.7) Therefore both eigenvalues of J are negative, indicating stable passive states for small Cr and confirming our numerical observations.The stability of passive states is lost when one or both eigenvalues of J become positive, therefore det(J) = 0 provides an equation to determine the critical Cr * .After simplification, we have Cr * k (1 − cos kπCr * ) = βr 2γ . (5.8) The root of Eq. (5.8) can be determined numerically, given that the parameters β and γ are known.In laboratory experiments and geophysical plate tectonics, the surface flow speed scale β and ventilation coefficient γ can be properly estimated and used to determine Cr * , therefore Eq. (5.8) offers the possibility of determining plate mobility through physical parameters.
We note that Cr * is usually small for large Γ as shown in Fig. 7, and the root of Eq. (5.8) in this limit can be shown as (5.9) Here we have used the relations r ∼ Γ −1 and k ∼ Γ , with the latter indicating the number of convection rolls is proportional to the aspect ratio Γ .Equation (5.9) can indeed be verified-as shown in Fig. 7(d), the critical Cr * measured from DNS data does follow the -2/3 power law with Γ .Much more can be investigated through the dynamics of Eqs.(5.3) and (5.4), and future experimental studies can certainly be used to further address the interaction between the plate and convective fluid below.

Discussion
In this work, we have numerically and theoretically explored the mechanical and thermal coupling between a moving plate and a convecting fluid beneath it.Inspired by present and past works (Gurnis 1988;Whitehead & Behn 2015;Mao et al. 2019;Mao 2021), 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 center 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 DNS is conducted at a parameter regime similar to laboratory experiments (Zhong & Zhang 2005, 2007a,b), 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 the flow has two counter-rotating large-scale circulations, whose strengths are modulated by the location of the plate (Zhong & Zhang 2005, 2007a,b).The plate is also bounded by the solid boundaries, so it stops moving once touching the wall.Modeling the plate dynamics in this case is not trivial, and advanced tools such as stochastic variational inequalities (Huang et al. 2018) can serve as a mathematical means of analyzing such plate-boundary interactions.A periodic domain of convecting fluid is necessary to experimentally verify 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 only investigate the dynamics of a single plate, our numerical method is capable of handling multiple plates providing that their interactions can be properly modeled.We are currently investigating such interactions, which has led to even more diverse and unpredictable dynamics.For example, if two small plates each have Cr < 1/3 but their combined size reaches Cr > 1/3, 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 analyzing the direct numerical simulation results there, we wish to further develop our stochastic model 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 fully capture the dynamics of continental drifts, they can hopefully still offer some fluid mechanical insights into the geophysics of Earth.
Nonlinear terms like u • ∇θ and u • ∇ω in Eqs.(6.10) and (6.11) are computed pseudospectrally with a simple and efficient anti-aliasing filter (Hou & Li 2007).With given initial and boundary data, (6.7) can be solved first to obtain θ n , which is inserted in f n so (6.6) can be solved next.Finally, (6.8) is solved with the known ω n .
After solving for the flow and temperature fields, the plate acceleration can be determined as a p,n = − Pr m P ∂ 2 ψ n ∂y 2 (x, 1)dx.(6.12) The plate velocity u p,n and plate location x p,n can then be computed through a 2 nd order Adam-Bashforth method, At Γ = 4, there are typically 256 Fourier modes in the x direction and 64 Chebyshev nodes in the y direction, with a time step size of ∆t = 10 −6 .These parameters are tested to yield resolved and accurate numerical solutions.

Figure 1 .
Figure 1.Schematics of the moving plate and convecting fluid.The fluid domain is bounded between y ∈ (0, 1) and periodic in x, and the floating plate of width d has center location xp.

Figure 2 .
Figure 2. Motion of a plate floating on top of a convective fluid.(a) Trajectories of plates of various sizes.The covering ratio is Cr = d/Γ and the time t0 marks the beginning of plate motion.(b) The total displacement of the plates, where dp is defined through ḋp = | ẋp|.(c) Average plate speed vp = ⟨| ẋp|⟩ becomes high when Cr > 0.33 and reaches a maximum at Cr = 0.56.In these simulations, Γ = 4, ρ = 4, Ra = 10 6 , and Pr = 7.9.

Figure 3 .
Figure 3. Dynamics of floating plate with Cr = 0.125 [(a)-(d)] and Cr = 0.5 [(e)-(h)].(a) A typical snapshot of the flow and temperature distributions under the plate with Cr = 0.125.(b) The y-averaged temperature and vertical flow velocity corresponding to (a), the shaded area indicates the position of the plate.(c) The plate displacement is a random process.(d) Plate velocity is a random variable with mean 0, whose PDF is a Gaussian distribution as indicated by the histogram.(e) At Cr = 0.5, typical flow and temperature distributions showing 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 nonzero mean.Supplementary movies of these simulations are included.

Figure 4 .
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 nonzero transnational 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 Fig. 2(a).

Figure 5 .
Figure 5. Equilibrium values of the plate velocity and phase angle.(a) Dimensionless plate velocity ûp is 0 for small Cr , but becomes translational as Cr increases.(b) The phase angle ϕ.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.
Fig. 2(a)  and has a transition between the noise-driven and linear motions at Cr ≈ 0.3.Thus, this simple model captures all the key features of the full numerical simulation.Without noise, the critical behavior of the dynamical system Eqs.(4.7) and (4.8) can be further analyzed.For small Cr , it is easy to see that Eqs.(4.7) and (4.8) have u p = 0, ϕ = 2πN (where N is an integer) as equilibria, which are stable and reflect the passive state of plate motion.Increasing Cr , new equilibria appear at Cr * = 1/3.For Cr > Cr * , it becomes possible to have a nonzero plate velocity u * p = (β/π)û * p , where û * p = − sin πCr Cr sin ϕ * (4.9) and the equilibrium ϕ * can be determined from cos ϕ * = 1 − (2Cr ) −1 cos πCr .(4.10)These states represent the translation of the plate.We note that Eqs.(4.9) and (4.10) are only functions of Cr , and 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 β/π.The possible values of û * p and ϕ * are shown in Fig.5.For Cr < Cr * = 1/3, u * p = 0 and ϕ * = 2πN are the only possible equilibria which reflect the passive nature of the small plate that is always attracted by the surface flow sink.For Cr > Cr * , new phases appear as ϕ * = (2N + 1)π ± arccos [(2Cr ) −1 − 1](cos πCr ) −1 , which are solutions of Eq. (4.10) and become stable for large Cr .They indicate that the larger plate tends to