Convection Cells With Accumulating Crust: Models of Continent and Mantle Evolution

Numerical calculations of two component cellular convection with a lighter component diffusing down from the top are a simple model of the convection of rocky planets with a crust. The crust either is mixed down, floats along the top as blobs (continents separated by ocean basins), or forms layers. The calculations are for very viscous fluid with density variations from temperature and the lighter component representing the crust. If variations are approximately the same size, the lighter component collects along the top into H‐shaped blobs/clusters that have a flat midsection with two lobes on each end. Elevations are calculated for a free upper surface and the shape resembles continent and ocean floor topography with level interiors and thickening (mountains) at the edges produced by sinking at the margins. Between each pair of clusters, thermal and concentration boundary layers resemble ocean basins with spreading centers. Convection is unsteady but introducing internal decay of the lighter concentration produces steady flow. Internal heating produces similar results along with periodic drifting and merging of blobs like some geological cycles. The fact that these features arise without phase changes or viscosity variation implies that blobs of continent‐like crust might be widely found on rocky planets.

Early work focused on motions caused by the presence of continents. Elder (1967) developed a theory and a laboratory demonstration showing the possibility that an insulating float of fixed size lying on top of cooling fluid might develop a spontaneous drift. Knopoff (1969) emphasized that the drift speed of plates is remarkably constant, which was in accord with Elder's results. This led to more controlled studies with theory and floating laboratory heaters (Howard et al., 1970) and modified theories with laboratory experiments (Knopoff et al., 1972;Whitehead, 1972). These produced motions of floats like the Wilson/supercontinent cycle showing back-and-forth steady speeds in a chamber rather than the changing speed of a passive tracer above convection cells. The end of the tank provides a reflection boundary condition, so a float reflection is equivalent to a supercontinent formation and breakup. Other laboratory studies of floating insulators on convection cells (Mao et al., 2019;Zhang & Libchaber, 2000;Zhong & Zhang, 2005, 2007 verify the self-generated movement. A monopole convection cell lying under an insulating float makes such a steady speed, rising at the trailing edge of the drifting float and sinking (with tilting) under the leading edge (Whitehead & Behn, 2015). Finally, layered convection with an internally heated top layer also develops growing oscillations that potentially lead to cycles (Busse, 1978;Rasenat et al., 2006).
Clearly, continents become modified through time, and questions abound concerning compositional changes of continent, water, and mantle over the past 3.6 billion years (Ko et al., 2022;Korenaga, 2021). The inclusion of thermochemical effects into numerical models has advanced greatly in the last decade. Examples are the lower mantle influence on the production of plume cycles (Li & Zhong, 2017), evidence of its importance near the core (Ko et al., 2022), and a recent numerical model of mantle convection in a spherical shell with accumulating continent-like clusters (Rozel et al., 2017). Although simplified models extend over the entire span of Earth's history, high-resolution models presently extend back approximately a billion years (Cao et al., 2021) and partially duplicate the Wilson/supercontinent cycle. This is less than 30% of Earth's age, so the detailed evolution of continent form and the oscillation like a Wilson/supercontinent cycle is being approached.
The overall picture is that convection in the mantle conveys the continents to many of the sinking regions. Lowman and Gable (1999) pointed out that heating in the mantle under continents occurred more from the absence of subduction than from insulation. Without any question, numerous advances in numerical programs over the past 40 years have made great contributions to the knowledge of continent-mantle convection interaction. It is becoming difficult to separate which components are needed for the predictive skill and which are simply adjusted to fit the present earth.
In contrast to many of the cited studies, the present model is purely predictive and does not try to include all the complexity of our earth. Results show the features of this model alone rather than a fit to present data. The model's balance is between convection and the form of the continent, including self-adjustment of its thickness and length. The model involves the convection of a fluid with a compositional variation making light fluid float on top of simple convection. It might be called an Occum's razor approach as it is one of the simplest of models for showing the evolution of a crust under the influence of mantle (Rayleigh-Benard) convection. A slowly diffusing composition variation (representing crust) in a host material (mantle) diffuses down from the top and represents crust as a simple material of lower density with no difference between continental and oceanic crust. Convection in the body of the fluid interacts with and deforms the crust to make structures near the top. Calculations span wide ranges of governing parameters over long evolution times so that the evolution of crustal form along with any drift or oscillation is resolved. This model answers the following two questions: 1. Is the form of the crust like or unlike the continents? 2. Can the crust split apart and drift or is it steady?

Methods
The numerical calculation consists of two-dimensional double diffusive convection in a rectangular chamber (the same geometry used by Elder, 1967through Lowman & Jarvis, 1999 to explore the simples' cases first) using the simple Matlab code given in the supplementary material. Equations step forward in time, driven by buoyancy from temperature and composition variations. All other physical properties are constant. The fluid has a very large and constant viscosity (infinite Prandtl number). The rectangular chamber has free-slip boundaries and zero heat flow and concentration diffusive flow out of the two sides. The bottom temperature is warmer than the top, producing well-known Rayleigh-Benard convection cells. In addition, compositional concentration C that increases density is advected and diffused with a smaller diffusivity than temperature. A lower value of C is set at the top boundary, and C diffuses to produce a lighter component (like fresher water at the surface of the ocean.) The result is a layer of lighter components that is swept to sinking locations. This layer is a hybrid ocean-continent crust. There is only one type of crust in this model and no difference between ocean and continental crust (unlike Earth, where the difference is significant).
The temperature (T), concentration (C), and Stokes flow equations (velocity ) are transformed into a dimensionless form suitable for numerical computation using D for depth scale, 2 / for the time scale ( is thermal diffusivity), κ/D for the velocity scale, and the external temperature and composition differences for T and C scales. The non-dimensional equations are The equations have some additional terms to well-known convection and diffusion equations. Uniform internal heat generation is quantified by internal energy generation E. With E = 0, convection is driven by a hot lower boundary, but with ≠ 0 the lower boundary is set to zero heat flow. In addition, heat production by the lighter component (perhaps representing radioactive heating in continent material) is given by the term in Equation (1) multiplied by ho. Finally, the term in Equation (2) multiplied by drives interior concentration toward the dense value of C = 1. The Stokes flow equation is valid in the limit of Prandtl number (Pr = ∕ ) infinite, where is kinematic viscosity. It is The dot signifies scalar multiplication, and the unit vector is in the downward direction of gravity. Flow is incompressible so a streamfunction defines the lateral (x) and vertical (z) velocities = − and = for a right-handed coordinate system (Whitehead & Behn, 2015, use a similar, but left-handed code) and the vorticity = − is related to it by The governing parameters are the Rayleigh number = Δ 3 ∕ that quantifies thermal buoyant forcing, the composition Rayleigh number = Δ 3 ∕ that quantifies compositional buoyant forcing, and the Lewis number = ∕ that quantifies the ratio of thermal to compositional diffusivities. A geometric dimensionless number is the aspect ratio of the chamber L/D where L is chamber length. The variable Δ is the temperature scale. For zero H, it is the excess of the bottom boundary over the top, or else it is the conductive temperature scale from H. The variable Δ is the concentration difference between the bottom and top boundary. In addition, g is the acceleration of gravity, is the coefficient of thermal expansion, is the coefficient of density change by composition, is the kinematic viscosity, and is the composition diffusivity.
In each numerical run, convection cells are first initiated with dense concentration C = 1 everywhere. At each time step, Equations (1) and (2) are numerically advanced, and then and are calculated, using Poisson equation solvers for the curl of Equations (3) and (4). Two-time step procedures were tested. The first was with central-difference time steps and the second was a forward-differencing time step. They produced similar results. The cycle is repeated many times. After convection becomes steady, a lower concentration (C = 0) is suddenly imposed along the top boundary of the fluid to initiate a layer of slowly accumulating lighter fluid. This layer thickens slowly, and ultimately large and small concentration mix within the fluid, but with ≠ 0 , the lighter component remains near the top of the fluid whether convection is present or not. Although explorations were conducted over a wide range of dimensionless number values, results shown here have dimensionless numbers in ranges that illustrate most clearly the competing effects of thermal convection cells with the lighter crust material.
For Earth, the difference between the surface elevation of the continent and the ocean floor is one of the most distinct aspects of its surface. The continent's surface (where we live) extends a few kilometers above the level of the ocean floor because the crust has a lower density than the mantle. Therefore, an estimate of surface elevation for this model is useful. Since the numerical grid has a flat top boundary, an equivalent elevation is estimated by calculating the pressure along the top to find a surface elevation as done in many numerical models (Ziethe & Benkhoff, 2010). Here, pressure is only a proxy for surface elevation. Unlike utilizing data for the earth that use the hydrostatic equation along with density and gravity to estimate surface uplift, this model assumes that Equation (3) is valid everywhere because the assumption is that the density of the host material is much greater than changes due to composition or temperature. Therefore, the lateral component of Equation (3) can be used with a calculated flow field to calculate surface pressure by integrating this laterally.
Let us imagine that this dimensionless pressure is equivalent to a surface elevation if the top surface were allowed to move vertically. The fluid has a density 0 and the surface pressure can represent an equivalent surface elevation compared to Earth using the following entirely synthetic calculation. The size of the elevation can be found using a relation between a dimensional surface elevation (proportional to surface pressure) and a dimensional overturning speed of a viscous overturning cell (aspect ratio of one assumed). The balances between elevation, pressure, and viscous flow are The only role of density in the model is to propel the cells with speed . We can use this to estimate an equivalent elevation for the earth. The viscosity value to be used is the mantle value such as = 10 19 m 2 /s based on inferences from the sinking speed of subducted lithosphere (Čížková et al., 2012). The speed is the approximate speed for plates = 3 × 10 −10 m/s, and d corresponds to the mantle depth d = 3 × 10 6 m (whole mantle convection). This = 3.3 km and this is a suggestive scale for elevation used in figures. It must be stressed that the elevation is modeled and not to be applied in detail to any natural application.
The first examples in Section 3.1 have a moderate value of Rayleigh number Ra = 10 4 with = ℎ = = 0 that drives a fully developed and slowly changing convection cell along with the stabilization by the lighter component for values of CRa ranging between 0 and 5 × 10 4 . Only two-dimensional flows are considered since convection rolls are stable in this range (Busse, 1967). The ratio of temperature diffusivity to component diffusivity is either Le = 100 or 10. The chamber is generally 4 times longer than the depth, but chambers 2, 8, 12, and 16 times longer than the depth have also been used to check the dependence of cellular wavelength on chamber length. Wavelength is similar to fluid depth in all cases. Therefore, cells in these studies have wavelengths close to fluid depth times pi just like cellular convection without a crust (Newell & Whitehead, 1969;Segel, 1969). Finally, examples are shown for greater values of Ra and CRa. They required both a smaller grid size and smaller time step thereby taking hours to days to complete. Only a few dozen were completed. The runs progress in time with the patterns of lighter component accumulations taking many different forms (see Video S1 of Supporting Information S1). Here we limit our attention to flows up to approximately t = 1 because C = 0 continues to accumulate and slowly accumulate at a strong top layer as, in Mercury for example (Ziethe & Benkhoff, 2010). The convection starts with uniform C = 1 and convection cells become steady by t = 0.2. Then, the lighter component (C = 0) is imposed thereafter along the top row of grid points. It diffuses down, forming a layer near the top. The lighter component is quickly swept laterally by the top of the existing convection cells and accumulates over the convergence regions above sinking thermals with thin boundary layers remaining over the diverging regions above rising locations. Three different things happen.  component ("continent") regions with thickening and surface elevations like "mountain ranges" along the margins. Each cluster has an H-shape with an interior of relatively constant thickness next to two sidelobes. Sinking occurs next to the margins and the localized convergence produces the lobes. The lobes represent mountains. The broadening of the clusters from t = 0.35 and 0.5 are maintained by small vertical upwelling hot motion rising near the bottom under the center of the elevated regions. Between each pair of accumulations, hot thermal plumes rise to the surface producing an equivalent to an ocean basin with thin thermal and compositional boundary layers with the middles elevated and diverging like spreading centers. The density structure resembles the lithosphere and ocean crust. The basin elevation is lower than the cluster elevation like "ocean basins next to continents".

Runs With
At t = 0.55 hot thermals from the deep convection region under the middle of the clusters rise and split the clusters apart. This resembles the breakup of a supercontinent. By t = 0.65 the new clusters (continents) have formed from the old pieces. Each half has drifted away from the centers of the previous clusters and collected at the previous ocean basin centers. Another earthlike feature is that at both t = 0.5 and 0.65 the temperature within each lighter component (continent) region is colder than the temperature elsewhere at the same depth, like the cold base of continents.
Video S1 of Supporting Information S1 shows the entire sequence up to t = 1. The structures in Figure 1a  To compare the elevation profiles in Figure 1a with Earth, the elevation of the seafloor and land surface around the equator from GeoMapApp is plotted in Figure 1b. Naturally, the earth's profiles are also affected by other 7 of 18 processes such as hotspots, differing spreading rates, and multiple plates on a sphere so the likeness breaks down in detail. Despite this, to the lowest order, the model at t = 0.5 produces the bimodal elevation of Earth from continent elevation and ocean depth. It also crudely captures the mid-ocean elevation of spreading centers, especially in the Atlantic, and the sidelobes that resemble mountain belts like the Andes.
In summary, the results are.
• The high surface elevations at t = 0.35, 0.5, and 0.65 are not only relatively flat like continents (tabular) but they also have local regions of higher elevation at their margins like mountain ranges. • The high surface elevations have cold roots just as continents do. • The surface elevation in the middle of the spreading center basins resembles ocean spreading center elevations. • The "continent" material migrates laterally. The movement has one complete cycle but then becomes irregular later.

Assorted Values of
Setting the three values = 1, 4, 10 in Equation (2) keeps the lighter component confined to the top regions and leads to more than 10 times longer duration of the clusters. Typical results are shown in Figures 2a-2f. The flow in 2a-c becomes completely steady in time (Video S2 of Supporting Information S1), so it is possible to study the structure of temperature, concentration, and streamfunction in detail. Figure 2c shows contours of vorticity and this highlights the role of the lighter component clusters. They produce vorticity that is of the opposite sign to the vorticity of the convection cells nearby. An intense region of shear lies next to the margins of the blobs leading to the thickening associated with the H-shape in the brown crust in Figure 1a at t = 0.5. This is true even with a chamber aspect ratio of 16 (Figure 2 d, e) and (Video S3 of Supporting Information S1). Figure 2f shows the average temperature gradient (heat flow) out of the surface over time, illustrating how the flow becomes almost completely steady. At t = 1.5 a cell near the end suddenly vanishes and the entire convection readjusts with a small change in heat flow.
The steadiness indicates a balance between the stabilizing lightness of the surface concentration with the underlying convection cells. It is probably aided by the fact that steady roll convection exists in this parameter range for convection without a crust (Busse, 1967).
Other runs have heat production added to the lighter component (ho = 40) as a simple model of continent crust heat production. The lighter component splits apart and merges cyclically with a neighbor over time (Video S4 of Supporting Information S1). It appears that the heat production for lightness near the surface produces some localized heating that destabilizes the cold continent roots and initiates the splitting, but a full investigation has not been conducted. This oscillation crudely resembles the Wilson/supercontinent cycle on Earth (Figure 2h).  0 ≤ ≤ 10 5 . The value C max = 0.5 is the chosen metric. It is safely used because the proportion is relatively insensitive to that value; see the inserted panel in Figure 3. Although unsteady flow makes scatter in some cases, the trend is clear. The proportion is large for CRa <10 4 signifying wide oceans and small or absent continents and it decreases with increasing CRa toward large continents, and small oceans.

Convection With Uniform Internal Heating
Convection with uniform internal heating has similar flows. It was calculated with zero heat flux imposed as a boundary condition along the sides and bottom with the product of E and Ra measuring intensity. A variety of steady and time-dependent flows was found for various values of CRa with HRa = 2, 4, and 8 × 10 4 and = 10 . Figure 4 shows the coverage and different modes of flow.
The triangles with an apex on the top (called "end" in the list in Figure 4) have relatively small CRa, and the lighter component has negligible buoyancy and minor influence on the flow. Flow occupies one long cell with time-dependent sinking on one end ( Figure 5, Video S5 of Supporting Information S1). The thermal boundary layer is unstable, and cold protrusions form at the top and are swept toward the left end by the surface flow near that end.
The circles in Figure 4 indicate steady flows ( Figure 6) that occur to the right at the largest values of CRa. The lighter component is swept to the sinking sites and retards the flow locally and produces a wide cold area of sinking. Upwelling occurs midway between the lighter component and produces surface elevation. The warmest isotherm has a plume-like structure even though it is not fed by a hot bottom boundary layer. The clusters of crust have a flat midsection with two lobes on each end. These also are models of continent crust like the runs in Figures 1 and 2. Therefore, the type of internal heating does not change the form of the floating material .  (Figure 7). CRa is roughly the same size as Ra. The lighter component is swept to the bottom and splits, but the plume goes alternately one way and then the other. Heat transfer along the top (Nu) changes during the cycle (Video S6 of Supporting Information S1).
The plus symbol in Figure 4 indicates flow that has time-dependent widening and contracting of the lighter component accumulations accompanied by faster and slower speed, and temperature and surface heat flow fluctuation (Video S7 of Supporting Information S1). The form of the temperature and lighter component accumulations fields are like those in Figure 6. Figure 8 shows two parts of the cycle.
The x-symbol in Figure 4 indicates a periodic flow that is somewhat like the Wilson/supercontinent cycle of the earth. The start and end of one complete cycle are shown in Figure 9 and many cycles are in Video S8 of Supporting Information S1. The lighter component splits and each half drifts apart only to merge with a neighbor and make a new blob of the lighter component. This is followed by a repeat of the cycle.

Ra up to 10 6
To determine whether these forms exist at larger Ra and CRa, calculations are done with Ra up to 2 × 10 6 and CRa slightly larger. They require finer grids and shorter time steps and at Ra > 10 5 , the value Le = 5 is used to save computer time. To establish the steady and time-dependent flows equivalent to Figures 6-9, we set = 10 3 . The flows have essentially the same form as lower Ra flows both for convection cells driven by a hot bottom temperature and cells driven by internal heating. The boundary layers for both the lighter component and temperature are concentrated in progressively smaller regions as expected. The flow is relatively steady but there are irregular

Discussion
The H-shaped clusters of crust that have a flat midsection with two lobes on each end as shown in some of the Figures are common in our results and found both with and without internal heating. The shape is linked to the two signs of vorticity at the margins (Figures 2c and 10) that are related to a new second harmonic in the convection cells (Figure 1a at t = 0.5). Additional computations about the two signs of vorticity with "kinematic cells" are in progress that show the same effect. These are not ready for this publication. In some cases, especially with the internal heating present, the rising protrusions near the bottom of the fluid under the lighter component regions are strong enough to reverse the sinking and then split the light accumulations apart (Figure 1a at t = 0.35, Figures 2g and 2h, Figure 9). In other cases (Figures 2a-2e), the bottom isotherm protrudes up under the lighter component but does not cause splitting. The videos show this effect dramatically.
Although this is a very simple model of convection with a crust, at least three dimensionless numbers are involved, and this might not be the simplest cellular convection model of continents and ocean basins. Two of them, Ra and CRa, are obviously needed since composition and temperature play important roles, and the third, Le would be very large, (>>10 2 ) in planets and Earth. Here, larger values of Le require large computation time, although so far, the difference between the results with Le = 5 and 100 seems minor. If this type of calculation is used in the future, a detailed study that illustrates the useful range of Le would be valuable. The fourth dimensionless number is the extinction coefficient for the lighter component with ≠ 0 . The dynamics are somewhat artificial, but it is useful to employ because it clearly produces a longer duration of the surface crust leading to structures and steady flow. The fifth number is the internal heating of the lighter component ( ≠ 0) and this produces cyclic flow. The sixth is uniform heating everywhere, E, and this also produces cyclic flow. Therefore, the number of relevant dimensionless numbers is closer to 4 or 5.
These calculations can be described to involve double diffusion and the new shapes found here can join the list of double-diffusive flow patterns like salt fingers and layered convection. Double diffusion possesses a great pedigree, having been discovered in the late 1950s from oceanographic questions about salinity variations in thermally stratified flows with more recent applications to magma cooling reacting flows, biology, and chemical processing (Huppert & Sparks, 1984;Radko, 2013;Schmitt, 1994;Turner, 1985). Despite the intense study of double diffusion, this application to mantle convection seems to be new since double diffusion is not even mentioned in three excellent textbooks (Davies, 2001;Schubert et al., 2001;Turcotte & Schubert, 2002). Therefore, the present results give both a simple new method and new perspectives about the dynamics of the distribution of surface crust with applications to planets, moons, continents, and ocean basins. A similar calculation using two immiscible fluids without double diffusion might be better as it removes the need for a specific value of both Le and Unfortunately, two-component numerical calculations are more challenging because the two materials can break up into small pieces or form very long thin strips and these add new numerical challenges as time progresses that are not present with diffusion. These present calculations also have the advantage of simplicity; the supplementary material contains the code in Matlab. The code with everything defined and integrated in time has only 118 important lines. Therefore, it might be difficult to find a much simpler Many studies of numerical mantle convection have included continents that are characterized by processes, not in this model. A table by Bobrov et al. (2022) lists 32 contributions and there are a small number of additional ones, too (i.e., Whitehead, 2017;Yoshida & Yoshizawa, 2021). Following the studies of drifting heaters and the suggestion that radioactive material in the continent contributes to drift (Elder, 1967through Whitehead et al., 2011, the next process hypothesized by many (including myself) to contribute to the Wilson cycle is that a rigid continent serves as a thermal insulator. The studies with insulating models of continents have been of gradually increasing complexity. Many of them show or analyze a partial cycle (Bobrov & Trubitsyn, 2008;Butler & Jarvis, 2004;Honda et al., 2000;Korenaga, 2007;Lenardic et al., 2011;Lowman & Gable, 1999;Lowman & Jarvis, 1993Rykov & Trubitsyn, 1996;Trubitsyn & Rykov, 1995;Yoshida, 2019;Zhang et al., 2009Zhang et al., , 2018 and some have achieved full cycles (Gurnis, 1988;Phillips & Bunge, 2007;Whitehead, 2017). All results indicate the universal nature of continent-convection interaction, and many quantify the role it plays in the evolution of Earth. A weakness of such models is generally that below the continents the temperature is elevated, which is well-known to be untrue. Likewise, models have come forth with continents of different viscosity or even with nonlinear rheology and some of these have also achieved a complete cycle (Bobrov & Baranov, 2019;Bobrov et al., 2022;Kameyama & Harada, 2017;Figure 7. Flip flop flow three times during a cycle. Ra = 10k, CRa = 12k, E = 8, = 10, ℎ = 0 . Video S6 of Supporting Information S1.

Figure 8.
Half-cycle of a periodic flow at the two times shown by the arrows. The widths and speeds increase and decrease. The total heat flow (red) is also cyclic (average heat flow = 4 but a value of 3 is added to the curve to get it on the same axis. (Ra = 10k, CRa = 15k, E = 8, = 10, ℎ = 0)( Video S7 of Supporting Information S1). Lenardic et al., 2000;Rolf et al., 2012;Yoshida, 2013;Yoshida & Santosh, 2011). The recent models examine the stress field during continent-to-continent collisions and breakups, including sea level changes of continent surfaces.
This present model is in contrast with all these increasingly developed models. It has no variation of thermal conductivity or viscosity to produce collisions, breakups, or a cold root. Results point to the roles of both seafloors spreading and subduction to continent material convergence. Recent models suggest that a large viscosity variation is a requirement to maintain the roots of cratons (Rolf et al., 2018;Yoshida, 2010aYoshida, , 2010bYoshida & Yoshizawa, 2021). This model cannot address issues of production, conservation, or redistribution of content material because the composition follows a diffusion/dissolution cycle unlike earth, so composition is not conserved.
Apart from Earth, the role of crust/convection interaction is important to study because the various rocky planets and moons have a range of crustal thicknesses. Only earth has a crust that deforms into continents. One might ask "is any resemblance between these model results and the earth, planets, or moons just a coincidence"? There are good reasons to think not. The diffusion of the lighter component here represents the differentiation of crustal material from the mantle (Spencer et al., 2017). If we were to build into any numerical code a dependence upon known differentiation processes (chemical reactions, liquid formation, porous flow, compaction of melt, and surface erosion), they would all spread out at different rates than the thermal field. These differences between the thermal time scale and the time of spreading of material are captured to the lowest order by our single Lewis number. Therefore, it is hoped that double diffusion calculations like these will become helpful to the understanding of mantle convection which includes crustal differentiation, continent formation, and interaction. Theoretical work still needs to be done. The three-dimensional counterparts to the two-dimensional flows given here are presently unknown, so extensions to three-dimensional cartesian bodies and spherical shells are encouraged.

Conclusions
Naturally, such a simple calculation has both strengths and weaknesses as a model of the mantle of earth and continents. The most interesting new effects seem to occur if Ra and CRa are similar in size, and they suggest how thick surface crust (continents) and mantle convection might interact.
The following strengths (with "best" parameters to make a form like the continents) are.  • Higher surface elevations occur over the lighter component at sinking regions like the continents. • The lighter component is wide with an interior plateau. • Some plateaus have higher elevations near the sides with deep roots like mountain belts near subduction locations. • The lighter component is colder than the convection at the same depth. This is like the cold roots of continents. • Sinking does not occur under the lighter component but is preferentially located next to them like Earth's continents. • The lighter component retards sinking and is responsible for a warm region near the bottom. • Early sinking of the lighter component reaches great depth. Figure 1a at t = 0.30 makes one ask, "are there implications for diamonds?" • The low elevation regions mimic ocean basins with a top cold thermal boundary layer and an oceanic crust that is advected laterally from the middle to the flanks of the basin. • Thermals eventually form under the "continents" and rise to split them apart (Figure 1a at t = 0.55). Does this resemble events when continents are split apart, and new spreading centers are initiated as in Africa today and at the beginning of the present Wilson cycle? • Internal heating added to the lighter component produces a cycle that is like the Wilson cycle. • The drift as evidenced by oscillation can be understood without viscosity variation. This is in accord with many of the early studies that had constant viscosity and produced drift easily. • The same drift/oscillation originates from deep heating under continents without continent insulation or continent heat production.
There are many important processes missing from this compared to full numerical models. A few of the following differences from duplication of Earth are • The lighter component blob occupies a large percentage of convection depth at Ra = 10 4 and this is unrealistically deep for our continents on Earth. Results for Ra = 10 6 are closer to the actual depths for Earth. • The spreading centers in Figure 1 come from deep plumes rising from the bottom boundary layer rather than only a divergence of surface mid-ocean rigid plates. Earth does not have sheetlike plumes under each mid-ocean ridge. This is not true with internal heating that still produces spreading centers. • Although there is cyclic splitting and merging of pairs of accumulations in some cases, there is no clustering of many to one supercontinents and therefore no complete Wilson cycle along with breaking apart and "continental drift" of many continents. • The constant viscosity everywhere is unrealistic. • The contribution of continent erosion is completely missing. It is suggested to fundamentally contribute to the continental crust thickness on earth (Whitehead & Clift, 2009;Whitehead, 2017;Wise, 1973 and references therein). • Earth clearly has heat flow from the core, internal heating, and possibly is still slowly cooling. Here, results have either a bottom hot temperature or internal heating but not both. • The ocean crust and continental crust are identical here, but they differ in composition on Earth.
It is vital to stress that this Occum's razor model is intended to show basic features of convection with a crust and not a detailed model of Earth; it is simply too reduced in complexity! Instead, only low-density material (a buoyant cluster) is included near the top of convection. It produces several features like Earth. Convection generates a wide tabular cluster with a cold root. The addition of internal heating produces drift but has the same cluster shape. The buoyant cluster has a cold base because the material in the cluster stays along the top for a longer time than the more rapidly moving fluid. Physically the cluster can be regarded as old as it is cooled by upward conduction over a long period of time compared to convective overturn. This makes the clusters resemble the continents. There is no need for continent insulation and viscosity variation. An important conclusion is that the vorticity produced by a cluster of light materials is important. It generates H-shaped clusters with lobes equivalent to mountain belts whether there are internal heat sources or not. It does this by producing localized stress near the cluster margins, reminding one of the deep-focus earthquakes that were the very first evidence that subduction zones contribute to mantle convection.

Conflict of Interest
The authors declare no conflicts of interest relevant to this study.

Data Availability Statement
All results are found using the MATLAB code like the one in Supporting Information S1. Figures 1 and 2 literally use this code. Figure 1b was obtained directly from GeoMapApp, which is available to the public. No additional data are included in this work. The supporting information also contains all eight videos.