Land Fraction Diversity on Earth-like Planets and Implications for their Habitability

A balanced ratio of ocean to land is believed to be essential for an Earth-like biosphere and one may conjecture that plate-tectonics planets should be similar in geological properties. After all, the volume of continental crust evolves towards an equilibrium between production and erosion. If the interior thermal states of Earth-sized exoplanets are similar to the Earth's, one might expect a similar equilibrium between continental production and erosion to establish and, hence, a similar land fraction. We will show that this conjecture is not likely to be true. Positive feedback associated with the coupled mantle water - continental crust cycle may rather lead to a manifold of three possible planets, depending on their early history: a land planet, an ocean planet and a balanced Earth-like planet. In addition, thermal blanketing of the interior by the continents enhances the sensitivity of continental growth to its history and, eventually, to initial conditions. Much of the blanketing effect is however compensated by mantle depletion in radioactive elements. A model of the long-term carbonate-silicate cycle shows the land and the ocean planet to differ by about 5 K in average surface temperature. A larger continental surface fraction results both in higher weathering rates and enhanced outgassing, partly compensating each other. Still, the land planet is expected to have a substantially dryer, colder and harsher climate possibly with extended cold deserts in comparison with the ocean planet and with the present-day Earth. Using a model of balancing water availability and nutrients from continental crust weathering, we find the bioproductivity and the biomass of both the land and ocean planet to be reduced by a third to half of Earth's. The biosphere on these planets might not be substantial enough to produce a supply of free oxygen.


Introduction
In the absence of observable unique biomarkers, the search for life beyond the solar system focuses on the search for habitable planets (e.g., Segura & Kaltenegger, 2010). The concept of the circumstellar habitable zone (CHZ) around a host star has long been established (Hart, 1979;Huggett, 1995) and developed further (see e.g., Cockell et al. 2016 andRamirez 2018 for reviews). The CHZ is defined as the zone of orbits in which a planet would have surface temperatures that allow the presence of liquid water, thus recognizing the importance of water for life as we know it (e.g., Brack et al., 2010). The boundaries of the habitable zone depend on the luminosity of the star but also on the atmosphere pressure at the surface of the planet and on atmospheric chemistry, in particular, on the concentration of the atmosphere greenhouse gases carbon-dioxide and methane (e.g., Lammer et al., 2009). More recent discussions of planetary habitability go beyond stellar luminosity and atmosphere chemistry to include other planetary properties such as the ratio between land and water covered surface areas, diversity of geologic provinces, magnetic fields, and processes such as plate tectonics that drive the planetary engine and its evolution (e.g., Heller & Armstrong, 2014;Schulze-Makuch et al., 2020;Glaser et al., 2020;Lingam & Loeb, 2018.
Having both emerged continents and oceans on Earth is an outcome of plate tectonics.
About 40% of its surface is covered with continental crust of which 87% is emerged; 13% of the area is continental shelf area and covered by water. Emerged land provides easy access to solar energy and weathering of continental crust rock provides essential nutrients for the biosphere. As a result, most of the biomass is found on continents and their shelves while the deep ocean biomass is much lower (e.g., Behrenfeld & Falkowski, 1997;Gray, 1997;Kallmeyer et al., 2012), even though the net primary production NPP from oceans is similar to that from land (Field et al., 1998).
The importance of a suitable balance between land and ocean covered surface areas for life as-we-know-it has recently been emphasized by Lingam & Loeb (2019 and Glaser et al. (2020). The required balance is a consequence of the need for water on the one hand and nutrients such as phosphorous and nitrogen on the other. On a land dominated planet with minimal water, the biological productivity should be low due to water limitations.
Ocean-dominated planets could be stymied by the lack of nutrients. If phosphorous is the productivity limiting nutrient (e.g., Kamerlin et al., 2013) then the increased weatherability of subaerial continental rock by slightly acidic rain-water gives land surfaces an edge over the ocean floors. The build up of free oxygen in the atmosphere for more advanced life is likewise favoured by continental crust weathering that produces the sediments in which organic carbon is buried and removed from oxidization. Lingam & Loeb (2019) discuss a simple parameterization of the NPP and the biomass as functions of the land-ocean surface ratio. They find both functions to be maximized at the present-day Earth ratio and to decrease by an order of magnitude as one goes to land dominated (1-2σ) or ocean dominated worlds. We note that the Earth-like land-ocean ratio would tend to maximise the extent of continental shelf areas where life will profit from solar insolation and the riverine entry of dissolved nutrients. Lingam & Loeb (2019 also find that oxygen can build up in worlds that have a water coverage between roughly 25 and 85%. Thus, planets with a balance of landmasses and oceans might be potentially conducive to hosting rich biospheres. In addition to the land-ocean distribution, the surface temperature will matter for the bioproductivity as has been emphasised recently by Heller & Armstrong (2014), Lingam & Loeb (2018) and Schulze-Makuch et al. (2020).
It is well established that the concentration of carbon-dioxide in the Earth's atmosphere is buffered by the long-term carbonate-silicate cycle (e.g., Kasting et al., 1993), thereby stabilising the atmosphere temperature and clement conditions. Negative feedback in the cycle causes the atmospheric CO 2 to attain an equilibrium between carbon outgassing and recycling (e.g., Lehmer et al., 2020;Isson et al., 2020). Plate tectonics is a central element of the modern cycle, exchanging CO 2 between the atmosphere and the mantle reservoirs through volcanism and subduction and the rate of outgassing is a function of the interior properties of the planet. The surface temperature is buffered as the rate of removal of CO 2 from the atmosphere through rain and weathering of continental crust -the latter a product of plate tectonics -increases with increasing temperature (e.g., Walker et al., 1981;Kasting et al., 1993;Foley, 2015). Weathering of continental crust results in carbon fixation and its rate depends on the continental surface area and the topography. Water is the main agent for transporting CO 2 from the atmosphere to the surface and for transporting weathering products including nutrients from the continents to the oceans and eventually to the subduction zones, where CO 2 is then returned to the mantle. Like CO 2 , water is cycled by plate tectonics through the mantle and mantle outgassing may have been a major source of the present ocean mass (e.g., Elkins-Tanton, 2011). The oceans provide the single most important thermal reservoir for buffering the surface temperature (e.g., Bindoff et al., 2007;Cheng et al., 2019;von Schuckmann et al., 2020) as well as they provide a sink for atmospheric CO 2 on short timescales (e.g., Pierrehumbert, 2010).
The operation of plate tectonics, like on present-day Earth, is vital to the long-term carbon cycle outlined above. Venus may have had surface recycling (e.g., Crameri, 2017) but its hypsometric curve is unimodal (e.g., Rosenblatt et al., 1994). While there are units resembling terrestrial rifts and ridges such as Beta Regio, there is no clear indication of subduction zones and continental crust, although Ishtar Terra as a high elevation physiographic province resembles a terrestrial continent (Ivanov & Head, 2011). Rather, it seems more likely that Venus has a continuous lithosphere. While Venus may be debated (Davaille et al., 2017), the smaller planets Mars and Mercury are widely agreed to have continuous outer shells termed stagnant-lid lithospheres. Planetary mantle convection modelling finds stagnant lids to form naturally as a consequence of a temperature dependent rheology of mantle rock (e.g., Stein & Hansen, 2008) while subduction and plate-like behaviour require additional strain related mechanisms for which water may be important (e.g., Bercovici, 2003).
Uncertainties in estimates of the interior properties of exoplanets have significant effects on assessments of their habitability (e.g., Seales & Lenardic, 2020). For example, uncertainties in the concentrations of heat producing elements will naturally affect estimates of the outgassing rates and the evolution of the surface temperature (e.g., O'Neill et al., 2020;Oosterloo et al., 2021;Unterborn et al., 2022). It has even been suggested that there may be feedback between the surface temperature and the tectonic regime (plate-tectonics versus stagnant lid) that could result in irreversible bifurcation between the two modes (Weller et al., 2015;Weller & Lenardic, 2018). Another key uncertainty in exoplanet evolution is the initial mantle temperature. Theoretical considerations and simple thermal history models suggest that differences in the initial mantle temperature should become insignificant after ∼1 Gyr of evolution due to the temperature-dependence of mantle viscosity (e.g., Tozer, 1967), at least for planets less massive or as massive as Earth (Kruijver et al., 2021). Therefore, one might conclude that the climate at 4.5 Gyr should be largely independent of the initial mantle temperature. This conjecture, however, neglects feedback associated with the growth of continental crust, which we will explore in this paper.
Even for Earth-sized planets with plate tectonics, it is not a given that the surface would have a balanced ocean to land coverage similar to Earth. Höning & Spohn (2016) and Höning et al. (2019a) have argued that -depending on initial conditions -the evolution of a plate tectonics Earth-like planet would rather result in a planet mostly covered by continents or mostly covered by oceans (Fig. 1). In their analysis, the continental crust and the water cycles are feedback cycles with either one or three fixed points depending on whether negative feedback or positive feedback dominates, respectively. Of the three fixed points, two would be stable (pertaining to the continent and the ocean planet) and one conditionally stable (the Earth). These differences in the outcome of planetary evolution are likely to be mirrored in the climate and possibly the biomass. The present climate of such a planet would retain a memory of the planet's very early evolution, since the equilibrium between CO 2 outgassing and silicate weathering would depend on the fraction of emerged The dependence of the strengths of positive and negative feedback on details of continental growth mechanisms complicates an extensive analysis of evolution pathways of (exo)planets (Höning & Spohn, 2016). Therefore, in the present paper we will parameterise feedback using four dimensionless quantities. We will show evolution pathways that depend on the assumed values of these quantities and then comprehensively discuss reasonable parameter choices and evolution scenarios. We extend our earlier models by analysing the effects of insulating continents delaying mantle cooling and the opposing effect of a depletion of mantle radioactive isotopes with continental growth. Furthermore, we calculate average surface temperatures for the model planets using a model of the carbonate-silicate cycle and estimate the bioproductivity and -mass following Lingam & Loeb (2021). We will present our models in a 3D phase space spanned by continental coverage, mantle water concentration and mantle temperature as variables.
In the following section, we begin by discussing plate tectonics from the point of view of complex non-linear system theory. We will then discuss fixed points of the system and their stability depending on dimensionless parameters that describe the strengths of positive and negative feedback. Coupling the model to a thermal evolution model of the mantle, we will show how fixed points evolve with time and steer the evolution of the planet and its habitability (Section 3). This section is followed by a discussion and conclusions.
2 Mantle thermo-chemical evolution and continental growth as a nonlinear system Plate tectonics and its role in geology has been described in a large number of publications (e.g., Sleep, 2015). In a nutshell, the Earth's surface is covered by seven major plates that move horizontally across the surface driven by solid-state convection in the mantle.
The plates consist of oceanic crust and lithosphere and continental crust and lithosphere.
Oceanic crust is produced by pressure release partial melting of mantle rock underneath mid-oceanic ridges and is returned to the mantle in subduction zones. The subduction zones are the places where modern plate tectonics produces continental crust by partially melting oceanic crust, subducted sediments, and mantle rock in the presence of subducted water reducing their solidus temperature (e.g., Campbell & Taylor, 1983).
Continents, despite their longevity in comparison with oceanic crust, are continuously weathered and eroded and the resulting sediments are transported to subduction zones by surface water flow. In addition, continental crust is eroded through friction between the continental keel and the subducting slab. The rates of subduction and crust production are closely linked to the rate of mantle flow, which is governed by the rheology of mantle rock.
On the scale of mantle flow, the rheology can be parameterised by an effective viscosity. The latter depends non-linearly on mantle temperature and on the concentration of water and -to a lesser extent -of carbon-dioxide in mantle rock. In addition, the mantle viscosity is a function of pressure, which should be considered for planets more massive than 2-3 Earth masses (e.g., Kruijver et al., 2021). Since we restrict our analysis to Earth-mass planets, we neglect this dependence here.
Sediments play important roles in subduction zone processes and in this paper, they play a decisive role: The erosion rate of the continents is proportional to their area and thus provides negative feedback limiting their growth. As carriers of water, however, sediments contribute to water transport (Hacker 2008;Deschamps et al. 2013)  While mid-oceanic ridges are the most prominent source regions of mantle water entering the oceans and the atmosphere, subduction zones are the major sink regions for surface water. Today, about 1.8 · 10 12 kg of water enter subduction zones each year mainly in hydrous minerals in sediments and hydrated basaltic crust (Jarrard, 2003). Significant, but difficult to constrain (e.g., Höning et al., 2014) quantities of the subducted water are expelled at shallow depth, the remaining water -partly stored in serpentinites and other hydrous minerals -reaches the source region of continental crust rock where hydrous minerals break down and release water that is consumed upon melting. The remainder of the water is subducted deeper into the mantle where it is eventually incorporated into nominally anhydrous minerals (Bell & Rossman, 1992). As a consequence, mantle viscosity is reduced (Hirth & Kohlstedt, 2003) and the rate of mantle convection flow increased. The water is transported in the convection flow towards mid-ocean ridges where part of it is outgassed into the atmosphere and oceans.
Feedback associated with the growth of continental crust can also originate from its thermal properties. The continental heat flow per unit area of about 70 mW/m 2 is known to be significantly smaller than the oceanic heat flow of about 105 mW/m 2 (e.g., Jaupart et al., 2015). Moreover, the basal heat flow per unit area from the mantle into stable continental crust is only about 15 mW/m 2 . These observations suggest that the continents act as thermal blankets of the mantle reducing the heat flow and the mantle cooling rate and diverting heat flow from sub-continental to sub-oceanic areas. The growth of continents diverting heat flow to oceanic areas can cause immediate positive feedback (Höning et al., 2019a). But even a reduction of the overall cooling rate introduces a sensitivity of the continental growth rate to its past history.
An opposing effect to thermal blanketing of continental crust is the redistribution of heat sources during crust growth, which may provide negative feedback. Radiogenic elements, because of their large ionic radii and their valence states, are enriched in melts producing continental crust. The concentration of heat sources in present-day continental crust is debated, however. Estimates are based on measured heat flow values and concentrations in mantle and crustal rock (e.g., Jaupart et al., 2016). Uncertainties are due to uncertainties in the ratio between surface heat flow and mantle radiogenic heat production (Urey ratio) and the density of heat sources in the lower mantle. But it is estimated that the modern continental crust contains between 25 and 40% of the total planetary inventory of heat sources (Dye, 2012;Jaupart et al., 2015Jaupart et al., , 2016. The transfer of heat sources from the mantle to the growing continental crust effectively cools the mantle (e.g., Spohn, 1991) and provides negative feedback by slowing down mantle convection and thus the production rate of new continental crust.
The above description of plate tectonics suggests that the continental crust volume, the mantle water concentration and mantle temperature are linked through non-linear feedback cycles as illustrated in Figure 2. An increase in mantle water concentration (blue) leads to a reduction of the mantle viscosity and thereby to an increase in the convection speed.
The increased subduction rate will increase the mantle water concentration establishing a positive feedback. The latter may be balanced by negative feedback associated with mantle water outgassing. The water cycle is linked to the continental crust cycle (green) through the water subduction rate. An increase in continental coverage increases the weathering and erosion rates and -assuming that the increased sedimentary volume will be subducted and contribute to water subduction -to increased rates of mantle water regassing and continental crust production. This feedback may be balanced by continental loss through erosion. The thermal evolution of the mantle (red) is linked to the mantle water cycle via the mantle viscosity and to the continental cycle via thermal blanketing and mantle depletion of radiogenic elements. The role of sediments in water subduction and the insulating effect of continents determine the strength of positive feedback in the coupled system.
The mantle temperature is subject to the negative feedback provided by the temperature dependence of the mantle viscosity and known as the Tozer principle (Tozer, 1967).
Accordingly, an increase in mantle temperature will decrease the viscosity which will result in an increased convection and mantle cooling rate -all else being equal. There is no positive feedback mechanism that could lead to runaway heating as long as the mantle heat source density decreases with time but continental crust growth and the mantle water budget can feedback onto the evolution of the mantle temperature. An increased degassing of the mantle will increase the mantle viscosity and result in reduced mantle cooling -with effects on crust growth and the water budget. Likewise, an increase in continental surface area can reduce the surface heat flow and the cooling rate -with possible readjustments of the convection system that may lead to adjustments of the plate speed. But only in an extreme case of thermal blanketing can continental growth result in temporary mantle heating, as we will show further below.
Mathematically, the non-linear feedback system can be described by the following set of differential equations. In writing Eqs. (1) and (2)   compare Eq. 1), mantle water concentration (blue, compare Eq. 2), and mantle temperature (red, compare Eq. 3). The evolution of the former two is coupled through the subduction of water. Positive feedback -indicated by '+' signs -results from the subduction of water reducing the mantle viscosity and increasing the plate speed and the crust production rate. Increasing the continental coverage results in an increased production rate of sediments feeding additional water and rock into the melt zone from which continental crust originates. Water outgassing, surface erosion and cooling dampen the positive feedback and are indicated by '-' signs. Continental crust growth feeds back further on the mantle temperature through thermal blanketing and transfer of radiogenic elements from the mantle to the crust. For further explanation see text.
-10-manuscript submitted to Astrobiologẏ rel. surface heat flow rel. mantle heat prod. +q core . (3) A, w and T are the state variables of the system. A is the continental surface fraction (combined surface area of the continents divided by the surface area of the planet), w is the mantle water concentration, and T the mantle temperature;Ȧ,ẇ, andṪ are their rates of change with time, respectively.Ȧ l,E andẇ l,E are the present-day rates of continental loss by erosion and mantle water loss by outgassing (in mass per time), respectively, and f r is the present-day ratio between loss and gain rates (f r = 1, for equilibrium between continental erosion and production and between mantle water gain and loss). We generally use the asterisk to denote variables scaled to present-day Earth values and the index E to denote the latter. v p (T, w) is the plate speed and a function of T and w. It is taken to be proportional to the mantle convection flow speed. L oo,oc (A) is the combined length of ocean-ocean subduction zones L oo (such as the Marianas) and ocean-continent subduction zones L oc (such as the subduction zone at the Peru-Chile trench) and is a function of A.
The length of subduction zones is calculated from the the total length of continental margins as a function of continental coverage (see Appendix B for details). ∆T u is the temperature difference across the upper thermal boundary layer of mantle convection, q m,E is the presentday heat flow from the mantle, Q A E is the radiogenic heat production rate as a function of time of the present-day depleted mantle (compare Eq. C1 in Section Appendix C). The depletion of mantle heat sources due to continental crust growth is modelled by the term 1 + ξ 4 (1 − A * ), where ξ 4 is the present-day share of the Earth's inventory of radiogenic elements that is in the continental crust. q core is the heat flow from the core. ρ m , c p , and V m are the mantle density, specific heat capacity, and volume, respectively. ξ 1 to ξ 4 are parameters describing the relative weight of terms in the balance equations (compare Tab. 1). As we will discuss in more detail below, ξ 1 and ξ 2 weigh the strengths of positive and negative feedback in continental growth, and f r shifts the fixed points in the phase plane.
Eq. 1 calculates the growth rate of continental crust area from a balance of continental crust production and loss through erosion. One part of the growth is directly proportional to the continental surface area. This part is a consequence of the effect of sediments on continental crust production and provides positive feedback as illustrated in Figure 2. Its relative weight is given by the dimensionless factor ξ 1 , which may take values between 0 and 1. A second contribution to the growth rate is proportional to the plate speed, which depends on the water concentration in the mantle and the mantle temperature, thereby coupling Eq. 1 to Eqs. 2 and 3. Its relative contribution is given by (1 − ξ 1 ). Loss by erosion is composed of a surface erosion part depending on the continental surface area and provides negative feedback to the system. The remaining part of continental erosion is related to the friction in subduction zones (subduction erosion) and depends on the total length of ocean-continent subduction zones. It should be noted that the latter part is a function of the continental surface area.
The term describing mantle water regassing (first term on the RHS of Eq. 2) mirrors the production term in Eq. 1, as both are related to the subduction of water. The regassing term altogether provides non-linear positive feedback to the rate of change of mantle water concentration through the dependence of the plate speed on the latter (compare Fig. 2).
Mantle water outgassing described by the second term on the RHS of Eq. 2 is directly proportional to w * , which contributes negative feedback to the change of the mantle water concentration. A step-by-step derivation of Eqs. 1 and 2 is given in Appendix A.
Eqs. 1 and 2 are coupled to the energy equation (Eq. 3) via the plate speed. The latter equation balances the surface heat flow with the heat production rate in the mantle by radiogenic decay and the heat flow from the core to calculate the rate of change of mantle internal energy. The surface heat flow is coupled to the water balance through the non-linear dependence of the mantle convection speed on the water concentration (Appendix C) and is calculated in terms of the present-day mantle heat flow −q m,E . Furthermore, it is coupled to the surface area of insulating continental crust through the term The term is written such that it will equal unity for A equal 0.4, the value for the present-day Earth, independent of the choice of ξ 3 .
ξ 3 is defined as the present-day ratio between the oceanic heat flow and the mantle heat flow. It is a measure of the efficiency of the thermal blanketing of the mantle by continental crust. We discuss a case of perfectly insulating continents, for which ξ 3 = 1,

Results
The state variables A, w and T span a phase space -a space in which all possible states of the system are represented by unique combinations of the state variables -in which solutions can be discussed. Equilibrium conditions within the phase space are described bẏ A = 0 andẇ = 0. These conditions define two surfaces in phase space that intersect in lines towards which the system will tend to evolve.
We begin by introducing the continent-water system in a phase plane evaluated at a representative present-day mantle temperature (Section 3.1). Following this, we discuss how the phase planes and the state variables evolve with temperature in 3D phase space in Section 3.2. In addition, we show time histories of selected quantities. Implications for the long-term carbon cycle are presented in Section 3.3 and the effect of thermal blanketing on the system is shown in Section 3.4.
3.1 Phase planes at constant present-day mantle temperature, fixed points, and stability and therefore shifts the fixed points in the phase plane. w and graphs of the integrals ofȦ = 0 (red) andẇ = 0 (blue). In panels a and b, with ξ 1 < ξ 2 these intersect in a single fixed point marking a single stable solution to the system of Eqs. 1-3 at constant mantle temperature T . The integral ofȦ = 0 is a single valued function A(w) under these conditions. Increasing the water content in the mantle will increase the growth rate through its dependence on v p (T, w). The increase is balanced by the erosion rate increasing with the continental surface area. The integral w(A) ofẇ = 0 is a single valued function for all parameter combinations considered with a maximum value around 50% continental coverage. The rate of water release first increases with mantle water concentration and is balanced by increasing water consumption in subduction zones. At around 50% continental surface area, the curve bends over as continents begin to dominate on the surface. For geometrical reasons, ocean-continent subduction zones and ocean-ocean subduction zones will shrink in length while the length of continent-continent collision zones without subduction will increase. This will reduce the potential to subduct water. The mantle water consumption rate will then decrease with increasing A and will need to be balanced, forẇ = 0, by a decreasing outgassing rate requiring decreasing values of w.
further discussion see the text.
with A for A > 0.37 is partly balanced by the sediment subduction flux increasing with A and the effect of the latter on mantle water regassing.
At ξ 1 ≥ ξ 2 , the integral A(w) ofȦ = 0 becomes a multivalued function of the mantle water concentration (compare Fig. 3, panels c and d). At small values of A, v p (T, w) increasing with w dominates crustal growth as before and A(w) has a positive slope on this first branch of the S-shaped curve. At a specific value of A, decreasing with increasing ξ 1 , the graph of A(w) bends over and its slope becomes negative. This occurs because the growth rate increasing with ξ 1 A can only be balanced by erosion if v p (T, w) and hence w decrease with increasing A. The increase of A(w) with negative slope continues until at some value of A, which is ∼ 0.6 in Fig. 3  The graphs of the integrals ofȦ = 0, A(w) andẇ = 0, w(A), intersect at three points (panel d) for ξ 1 > ξ 2 . Of these, two represent stable equilibria at small and large continental coverage, respectively. A(w) has positive slopes at these points. The intersection point between the two -at A = 0.4 and w = 100ppm here -is a saddle point, stable with respect to perturbations in w but unstable with respect to perturbations in A. A(w) has a negative slope at this point. In panel c, we plot a transitory case for which ξ 1 = ξ 2 = 0.5, f r =1. Here, the two curves just touch at A = 0.4 and w = 100ppm with dA/dw = dw/dA and there is a stable equilibrium point at A = 0.25 and w = 80ppm. It should be noted that ξ 1 = ξ 2 = 0.5 marks a bifurcation for the solutions of Eqs. 1 -2 (Höning et al., 2019a).

Mantle cooling without thermal blanketing
In the following, we will present numerical solutions to Eqs. 1 -3 in the A, w, T -phase space (Fig. 4). We also show phase planes at the present day mantle temperature and at 1850 K, respectively, with graphs of the integrals ofȦ = 0 andẇ = 0. The latter phase planes can be compared with those in Fig  Model calculations are started at t = 0 by integrating Eqs. 2 and 3 with a given initial temperature T 0 and water content w 0 = 100 ppm. At a later time t = t p according to the assumed onset of plate tectonics -or more specifically, the onset of continental growth by subduction zone melting -we include Eq. 1 in the integration with A 0 = 0 as starting value. The onset time t p is 0.5 Gyr for the models in Figs. 4 and 5. Note that the onset of continental growth by subduction zone melting needs to be started that early as we grow the entire continental surface area by that mechanism. The starting values could be easily adjusted if plate tectonics is taken to start with the Proterozoic at t p = 2.5 Gyr and the Archean continents had grown by a different mechanism to some area A 0 > 0 at the time.
We present three models in each panel of Fig. 4 that differ by the chosen initial mantle temperature, colour coded cyan, green and purple. The initial mantle temperature for the models belonging to the green curves (1945 K) in both panels was chosen such that a continental coverage of 40% as on present-day Earth results after 4.5 Gyr. For the models belonging to the purple curves, a 150 K higher initial temperature was chosen; for the models belonging to the cyan curves the initial temperature was chosen 150 K lower. Otherwise the models have the same parameter values. We note that at given values of ξ 1 , ξ 2 , and f r , the values of T 0 and t p can be traded against each other to arrive at the Earth's present-day continental coverage. For instance, we could have used a lower initial mantle temperature in combination with a later onset time t p . We could also have set the value of f r > 1, i.e. forced the present-day Earth continental growth rate to exceed the erosion rate, which would, for example, have allowed a smaller value of T 0 . conditions differences become small after about 2 Gyr of evolution. Again, the models with ξ 1 > ξ 2 evolve less rapidly with a similar delay as for the mantle temperature.
Differences between models with dominating positive versus negative feedback, respectively, are more significant for the mantle water concentration and, even more so, for the surface coverage by continental crust. With negative feedback dominating, ξ 1 < ξ 2 , the mantle water concentration adjusts rapidly to initial conditions and increases steadily after a little more than 1 Gyr almost independent of the choice of initial mantle temperature to reach a present day value of about 70 ppm, the equilibrium value being about 65 ppm. For positive feedback dominating, ξ 1 > ξ 2 , the initial adjustment of the variables to their initial conditions at times smaller than 0.5 Gyr corresponds to the one in the model with dominating negative feedback. But beyond 0.5 Gyr, the evolution paths diverge with greater and widening differences in mantle water concentrations and continental coverage. This is caused by the emergence of the three fixed points and the evolution paths tending towards them.
Figs. 4 and 5 illustrate how positive feedback causes the present-day value of continental coverage to depend on its early evolution. Although mantle temperature values converge after a few billion years of evolution and differences in initial conditions disappear as predicted by the Tozer principle (Tozer, 1967) imply a wet mantle while a small coverage is related to a comparatively dry mantle. This is due to the fact that both gain rates -continental production and mantle water regassingdepend on the same physical parameters that control the transport of water in subduction zones. Because we neglect the loss of water due to atmospheric escape, the complement of the mantle water is to be found at the surface as ocean water.

Implications for the long-term carbon cycle
The surface area of continental crust emerged above sea level has a direct impact on the climate because the rates of erosion and silicate weathering increase with increasing land area (e.g., Foley, 2015;Graham & Pierrehumbert, 2020;Levenson, 2021). A planet with a larger area of emerged continents will thus have a lower concentration of CO 2 in its atmosphere (all other factors being the same) and a lower surface temperature as compared with a planet with a smaller land area. For details of the calculation of the CO 2 degassing rate and partial pressure and the surface temperature see Appendix E. In Fig. 6, we plot the

Cooling with thermal blanketing and mantle depletion in radiogenic elements
Thermal blanketing of the mantle by the continental crust and the depletion of the mantle in radiogenic heat sources transferred to the crust can provide additional feedback between the crust growth rate and its surface area. In Fig. 7, we plot results for models with the same parameter values as in Fig. 4 but account for thermal blanketing by setting ξ 3 = 0.85 and ξ 4 = 1/4 in Eq. 3. The onset time t p of continental growth by subduction zone melting for the green model was adjusted such that the continental surface fraction grows to 0.4 in 4.5 Gyr. For models with ξ 1 = 1/3 and ξ 2 = 2/3, t p was adjusted to 0.253 Gyr. For those with ξ 1 = 2/3 and ξ 2 = 1/3, t p was chosen to be 0.161 Gyr.  Interestingly, for the case of dominating negative feedback (Fig. 8 left), large continents (purple curve) result in comparatively hot and dry mantles, with opposing effects on viscosity and plate speed. For dominating positive feedback (Fig. 8 right), sediments from the comparatively large continents contribute to water subduction such that a hot and wet mantle results, both increasing the plate speed. This feeds back to continental growth such that the difference in continental coverage at 4.5 Gyr becomes particularly pronounced. Figure 9 shows results pertaining to the long-term carbon cycle for the models that include thermally insulating continents. As in previous figures, panels on the left are for models with predominating negative feedback and those on the right are for models with predominating positive feedback. A larger surface fraction of insulating continents keeps the mantle warmer and therefore results in somewhat larger degassing rates but also in a larger weatherability and rate of CO 2 removal. The differences in the CO 2 partial pressure and surface temperature between the models remain small and are even a bit smaller than for models neglecting thermal blanketing altogether.

Discussion
In this paper, we modelled three main geologic processes as a nonlinear coupled dynamical system: growth of continental crust, exchange of water between the reservoirs on and above the surface (oceans, atmosphere) and in the mantle, and cooling by mantle convection.
These processes are linked through mantle convection and plate tectonics with • subduction zone related melting and volcanism and continental erosion governing the growth of the continents (Eq. 1), • mantle water degassing through volcanism and regassing through subduction governing the water budget (Eq. 2), and • heat transfer through mantle convection governing the thermal evolution (Eq. 3) The rates of these processes depend non-linearly on mantle temperature and water concentration through the mantle viscosity (compare Eq. C6) and plate speed v p (T, w) which couples the governing equations 1 -3. By including a model of the carbonate-silicate cycle, we explored consequences for the atmosphere's CO 2 partial pressure and average surface temperature as indicator of the planet's habitability. In addition to the coupling through v p (T, w), there are immediate feedback mechanisms with the crust production and erosion rates depending on the continental surface area and the length of continental margins and with the mantle water outgassing rate depending on the mantle water concentration. The former is implied by the dependence of the melt production rate in the subduction zone on the water supply through sediments (Höning & Spohn, 2016) which in turn is directly proportional to the continental surface area. The latter is implied by the dependence of the mantle solidus temperature and the viscosity on the water concentration. The mantle water concentration and the continental surface area feed back on the mantle thermal energy balance; the former through the dependence of v p (T, w) on the mantle water concentration, the latter through the thermal blanketing effect of the continents.
The conceptual nature of our study motivated us to simplify our model beyond the level of abstraction of our previous work (e.g., Höning & Spohn, 2016;Höning et al., 2019a) by reducing the number of free parameters to four, ξ 1 , .., ξ 4 , with 0 ≤ ξ 1 , .., ξ 4 ≤ 1. We still use the present-day Earth as a reference and present-day observables as scales. In the following, we will discuss our finding that the spread of continental coverage on Earth-like planets is determined by the respective strengths of positive and negative feedback in continental growth (ξ 1 vs. ξ 2 , Section 4.1) and by the relationship between thermal blanketing and depletion of radioactive isotopes upon growth of the continental crust (ξ 3 vs. ξ 4 , Section 4.2). Uncertainty in these parameter values represents the main uncertainty in the model.
Our discussion aims to provide a better qualitative understanding of the feedback processes; we admit to lacking data for a detailed understanding of quantitative differences. Finally, we discuss the habitability of the model planets and their bioproductivity that result from the various possible evolution pathways in Section 4.3.

The feedback parameters ξ 1 and ξ 2 and their implications for continental growth
The parameters ξ 1 and ξ 2 can be regarded as slide rulers in Eqs. 1 and 2 to weigh between feedback terms. We have chosen combinations of ξ 1 = 1/3 and ξ 2 = 2/3, ξ 1 = ξ 2 = 0.5, and ξ 1 = 2/3 and ξ 2 = 1/3 to demonstrate in Figs. 3 -9 how the properties of the system vary with these choices. Note that ξ 1 and ξ 2 do not have to add up to 1, although they do in our parameter choices. The parameter ξ 1 weighs between the importance of the continental surface area and plate speed (that is, convection speed), respectively, for both the production rate of continental crust and the mantle water degassing rate of the present-day Earth.
The value of ξ 1 becomes particularly important late in the evolution, when the mantle has cooled down sufficiently such that v * p , the plate speed scaled to its present-day value, will be of the order of unity. Geologically, this would be in the post-Archean. At earlier times, v * p will have a greater effect than A, unless ξ 1 is chosen such that (1 − ξ 1 ) 1. As a geological interpretation, ξ 1 weighs the importance of sediments for transporting water in the subducting slab to the source region of continental crust for the present-day Earth.
The rate of melting and crust production in the model is directly It may be argued that positive feedback may also originate from the variation of v * p with viscosity. After all, an increase in v * p will reduce the viscosity and increase the convection speed and, in turn, the water regassing rate. This will be counteracted, however, by increased water outgassing and cooling associated with the convective mantle flow. In the present model, the effect of cooling on increasing the viscosity dominates over the effect of an increasing water concentration. Altogether, the contribution of convection and plate circulation is to provide negative feedback. Other positive feedback mechanisms may exist that are independent of erosion and sediment subduction. We will discuss the contribution of thermal blanketing to positive feedback in the next section.
The parameter ξ 2 weighs between surface erosion and subduction erosion, the latter depending on the length of the ocean-continent subduction zones. A major share (≈ 1/3) of the global continental erosion rate has been associated with subduction erosion (Stern, 2011).
In addition, significant shares of the loss of continental surface area have been attributed to continental collision and lower crust delamination, leaving less than a third for surface erosion Stern (2011). It thus seems reasonable to assume that ξ 2 < 0.5. We note that the rate of subduction erosion may depend on the rate of friction between the subducting slab and the continental keel in addition to its dependence on the length of the subduction zone. Sediments may reduce friction by lubricating the slab Stern (2011) Our results show how strong negative feedback would lead to an evolution largely independent of the starting conditions and the early history of the planet which would imply a single stable present-day value of the continental surface area. For strong positive feedback, however, the outcome of the evolution may be quite different depending on starting conditions and the early history. For exoplanets this suggests that Earth-like planets covered largely by land (i.e. land planets) are possible as well as planets covered mostly by oceans (i.e. ocean planets) (see Fig. 1). Höning & Spohn (2016) found that the land planet should be the most likely outcome, followed by cases that end as ocean planets. Only a small percentage end as planets with a balanced Earth-like share of land and ocean surfaces.
Earth's further evolution can be discussed using the present day phase planes shown in Figs. 4 and 7. It will depend on the location of the state variable A relative to the fixed point.
For values of A close to but larger than the continental coverage at the unstable fixed point, A will increase with time at a small rate beyond its present-day value. While A is growing, the mantle water concentration will increase. With water budget conserved, this will be accompanied by decreasing surface water reservoir and increasing area of emerged land.
This would be beneficial for the biosphere allowing for more silicate weathering, provision of nutrients, and partly compensating for an increasing solar luminosity. Eventually, the unstable fixed point will overtake the trajectory and net continental growth will cease.
Additional calculations have shown that the upper stable fixed point will move towards smaller values of A together with the unstable fixed point moving towards larger values of A. The two fixed points will eventually merge and then disappear with further cooling while the continental area will slowly decrease with time being far away from but evolving towards the remaining lower stable fixed point.
An extreme case in which continental production and erosion would be entirely independent of the surface area can be studied by setting ξ 1 = 0 and ξ 2 = 0. The continent production rate will then depend only on the mantle convection rate and the erosion only on the total length of ocean-continent subduction zones. In addition, we neglect thermal blanketing (as in Fig. 4)  It is interesting to note that for planets around M-dwarfs, Tian & Ida (2015) predict a similar bimodal distribution of emerged land area, with most planets either having their surface entirely covered with water or with significantly less surface water than Earth. Their analysis, however, was based on a combination of planet formation and water escape models, not explicitly modelling continental growth.
The small positive net growth rate of the continents since the Archean (e.g., Armstrong, 1991;Dhuime et al., 2017;Rosas & Korenaga, 2018;Korenaga, 2021) provides a strong argument for positive feedback in Earth's continental growth history: Positive feedback can sustain a small net growth rate over an extended period of time even in the presence of a slowly decreasing mantle temperature as in the right columns of Figs. 5 and 8. In that case, the unstable fixed point will be evolving slowly towards larger continental coverage and A, being located slightly above this fixed point, will follow suit. We note that equilibrium is not expected to be attained as mantle heat production decreases with time and planets necessarily cool as has already been discussed by Schubert et al. (1979). With negative feedback prevailing, the continental production rate will decrease with decreasing mantle temperature. As a result, the continental surface area will decrease with time as in the left columns of Figs The continents may have grown in the Archean by mechanisms that differ from those of modern plate tectonics (e.g., Condie & Kröner, 2013). Our analysis would still be applicable, because the bifurcation appears only at post-Archean mantle temperatures and plate speeds.
The onset time t p could be 2.5 Gyr and a value for A at the time would be >0 along with assumed values for T and w. In such a scenario, the present-day continental coverage would be sensitive to the coverage at t p .
4.2 Thermal blanketing and mantle depletion through crust growth; effects of ξ 3 and ξ 4 Thermal blanketing by continental crust can have two effects that may both act at the same time: Blanketing can result in an overall reduction of the surface heat flow and the cooling rate. And it can result in a redirection of mantle heat flow away from the continents to the oceanic areas. Using stability analysis, Höning et al. (2019a) have shown that thermal blanketing can lead to runaway continental growth and even mantle heating (see also Lenardic et al. 2005 for related findings). The analysis, however, was based on the assumption of a steady state mantle and a constant Urey Ratio. It also did not include the redistribution of radiogenic elements between the mantle and the crust. In their analysis, blanketing caused the heat flow from the mantle to be diverted from the sub-continental to the sub-oceanic areas, resulting in an increasing sub-oceanic mantle temperature, heat flow, plate speed and, finally, continental growth rate. Their result can be considered as an extreme case but suggests that positive feedback between thermal blanketing and continental growth is, in principle, possible.
In the present model, we consider thermal blanketing in the presence of mantle cooling and the decay of radiogenic elements. In addition, we consider the transfer of heat producing elements from the mantle to the crust as a consequence of melting and crust growth. The cooling and the decrease in time of the mantle heating rate will reduce the self-reinforcement of continental growth through thermal blanketing. We do not consider any possible redirection of mantle flow away from sub-continental areas to sub-oceanic areas. Instead, we use ξ 3 as a slide ruler to distribute the present-day mantle heat flow q s,E between the oceanic areas and the continents. From its definition in Eq. D3, ξ 3 is the present-day share of the oceanic heat flow in the total heat flow from the mantle. For simplicity, this ratio is kept constant throughout the evolution calculation. Also, we neglect any enrichment of the oceanic crust in radiogenic elements. gives 11 -12.3 TW. However, this value is much too large to be explained by continental heat production for which Jaupart et al. (2015Jaupart et al. ( , 2016 give 6 -7 TW, the remainder of 3 -4 TW to be attributed to continental margins where the heat flow can be significantly higher than in stable continental areas (e.g., Springer & Förster, 1998). Altogether, this suggests a global present-day heat flow from the mantle of ≈ 38 TW, so that values for ξ 3 between 0.8 and 0.85 should be reasonably characterising the heat flow ratios on the present-day Earth.
To calculate a consistent value of ξ 4 , the share of the heat sources in the present-day continental crust (compare Eq. C3), we use the observed continental heat flow and subtract the above estimate of the mantle heat flow through the continents including continental margins of 6 -7 TW to obtain 7 -8 TW. Unfortunately, the mantle heat production rate is much more uncertain. Dye (2012) cites estimates between 11 and 24 TW for the bulk silicate Earth (mantle + crust) as resulting from cosmochemical studies. Using a mantle heat flow estimate of 36 TW and estimates from convection modelling of 30% to 50% being attributable to radiogenic heating (Dye 2012, the so called Urey ratio), they arrive at 19 to 27 TW, again for the bulk silicate Earth. Note that we do not adopt the case of 80% radiogenic heating as we find that the Urey ratio predicted from constant viscosity mantle models is unrealistic. Jaupart et al. (2015) have compiled mantle and crust heat production rates from the literature and from their own interpretations of global heat flow data and arrive at a representative value of 11 ± 6 TW for the mantle heating rate. This value can be considered as a lower limit as it assumes that the (depleted) MORB source extends through the whole mantle. Thus, values for ξ 4 between 0.25 and 0.4 should be reasonable to assume.
In Fig. 7, we used a combination of values for ξ 3 and ξ 4 within our range of estimates that maximises the effect of thermal blanketing. For values of ξ 4 > 0.25, depletion of heat sources works increasingly against the effect of thermal blanketing such that for ξ 3 = 0.8 and ξ 4 = 0.4, the results are close to those neglecting blanketing as in Figs. 4 and 5. In Fig. 11, we show the difference of the continental surface fraction at 4.5 Gyr resulting from models with a difference in the initial mantle temperature of 300 K (i.e. between the light blue and purple curves in Figs. 4 and 7) for dominating negative (left) and positive (right) feedback.
As expected, the spread in continental surface area at 4.5 Gyr decreases with increasing ξ 4 but the spread can be twice as large for cases with dominating negative feedback and ξ 4 = 0.25.
Long-term heating of the mantle can only be obtained if the thermal blanketing is made extreme. In Fig. 12, we show a case with ξ 3 = 1 implying perfectly insulating continents.
Here, with dominating positive feedback for a high initial mantle temperature, we obtain heating even though the depletion of radiogenic elements in continental crust is taken into account with ξ 4 =0.25.
All in all, we conclude that the reduction of cooling by thermal blanketing underneath growing continental crust is partly, potentially largely compensated by the effective cooling of the mantle associated with the transfer of radiogenic elements from the mantle to the crust.

CO 2 outgassing, weathering, surface temperature and bioproductivity
Coupling the continental erosion rate and the plate speed to a model of the long-term carbon cycle Appendix E, we find that the large spread in continental surface coverage between our models results in differences of surface temperatures of not more than about 5 K. The reason for this rather small difference are opposing effects of increasing continental coverage on the atmospheric CO 2 concentration: On the one hand, the weatherability and thereby the CO 2 sink increases with increasing continental coverage. On the other hand, the mantle remains hotter with larger, insulating continents, which causes an enhanced outgassing rate of CO 2 . Note that our ocean planet is a shallow water world with an ocean depth of about 4 km for which the carbonate silicate cycle is still operative (e.g., Höning et al., 2019b;Hayworth & Foley, 2020), even though the negative feedback associated with silicate weathering would be weaker. It should not be confused with a deep ocean water world with 10 -1000 Earth ocean masses for which the carbonate-silicate cycle might not work at all due to the formation of high-pressure ice layers and potentially limited magmatic degassing (e.g., Noack et al., 2016;Kite & Ford, 2018;Krissansen-Totton et al., 2021).
A 5 K surface temperature difference is considerably smaller than the excursions of the Earth's mean surface temperature over the Phanerozoic (e.g., Scotese, 2016) and might be assumed to be of limited significance for the general discussion of exoplanet habitability.
However, Lingam & Loeb (2018) and Schulze-Makuch et al. (2020) have argued that a 5 K higher surface temperature would significantly increase the habitability of an Earth-like world. In combination with the differences in land and water surfaces it could have an impact on climatic zones. Earth climate modelling suggests that the temperature difference on land surfaces, in particular near the equator, may be twice as large (Nazarenko et al., 2022).
Furthermore, the polar ice sheets will shrink with increasing temperature. In contrast, the lower surface temperature of the land planet would allow for substantial polar ice caps.
These differences should be reflected in the sea-level, which would also differ on account of the non-linear thermal expansion of seawater (e.g., Widlansky et al., 2020). Undoubtedly, the biosphere on worlds with different land fractions would adapt to their environment, and a balanced ratio of ocean-to-land accompanied by large continental shelves is certainly beneficial for an Earth-like biosphere (e.g., Kallmeyer et al., 2012;Höning & Spohn, 2016;Lingam & Loeb, 2019;Stevenson & Wallace, 2021). Vice versa, life on Earth co-evolved with the environment, shaping the atmosphere and oceans (Grenfell et al., 2010;Lenton et al., 2014;Algeo et al., 2016). The emergence of pelagic calcifiers in the deep oceans has been argued to have played a role in the stability of Earth's climate (Ridgwell et al., 2003) and the colonisation of land by plants enhancing continental weathering has presumably been crucial in keeping our planet's biosphere flourishing in the long term (Berner, 1998;Schwartzman, 2017). In combination, life on land and in the oceans has been shown to help stabilising Earth's climate on timescales >100 kyr with land plants dominating on timescales >1 Myr and marine calcifiers being predominantly important on shorter timescales (Höning, 2020). This particular relative importance of life on land and in the oceans is a direct result of Earth's balanced ratio of oceans to land, however: On planets with a larger continental coverage, land plants would become increasingly important for climate stability, whereas the role of marine calcifiers would become minor (Höning, 2020).

Conclusions
In this paper, we assessed the coupled evolution of continental growth, mantle water cycling, and thermal evolution of Earth-like planets. In particular, we explored how differences in initial mantle temperature and the strength of positive and negative feedback would impact their continental growth curves and water budgets. We coupled the model to a model of the carbonate-silicate cycle via the plate speed and the continental surface area to assess the habitability and the bioproductivity of these planets. Our results can be summarised as follows: • Positive feedback in the coupled continent-water cycle may lead to the emergence of multiple fixed points and to substantially differing land and ocean coverage, with an Earth-like distribution of land and ocean surfaces being the least likely. The average depth of the oceans would be similar, around 4 km, because differences in the ocean surface will be compensated by differences in the distribution of water between the mantle and the surface reservoirs. The final continental coverage and the environment of the planet will depend on the initial conditions, in particular on the initial mantle temperature and the onset time of plate tectonics. If early mechanisms of continental growth differ from plate tectonics then other differences in initial conditions may still matter. An evolution largely independent of the early conditions is possible only if negative feedback mechanisms prevail.
• Positive feedback is mostly provided by the role of sediments in continental growth.
The rate of sediment subduction will increase with the continental surface area and the additional water transported will enhance the continental growth rate. In addition, sediments lubricate the subducting plate and facilitate efficient subduction.
Thermal blanketing of the mantle by continents may retard mantle cooling, keep an enhanced continental production rate and may cause an even larger difference in final continental coverage. But the effects of thermal blanketing is partly, potentially largely, balanced by the removal of heat sources from the mantle going along with crust growth and with an enrichment of radiogenic elements in the continental crust.
• The carbonate-silicate cycle is coupled to the evolution of the interior through the carbon outgassing rate depending on the plate speed and to continental growth via the weathering rate increasing with the continental surface area. We found differences in the mean global surface temperature between the land planet and the ocean planet of up to ∼ 5 K. The climate of all three planets as well as their biospheres may have evolved quite differently, however. The warmer ocean planet's potentially higher specific bioactivity may be limited by the smaller availability of critical nutrients from weathering its smaller land surface. The cooler land planet may have a smaller specific bioactivity further reduced by the limited availability of water but its bioactivity would be spread over a larger area. By using the scalings of Lingam & Loeb (2018 we find the bioproductivity and the potential biomass on these planets reduced by a third to half of Earth's. It is questionable whether a biosphere on these planets would be substantial enough to produce more oxygen then is consumed by geologic activity.
• Main model uncertainties are related to quantifying the role of sediments in continental production and to estimating the depletion of radiogenic elements by continental crust production. Future work on these issues is required in order to improve understanding the complexities of our own planet and the prediction of the continental fraction on Earth-like exoplanets.

Appendix A Coupled continental and mantle water cycles
Key variables in the model are the surface fraction of continental crust A and its rate of changeȦ, which depends on the rates of continental productionȦ g (index g for gain) and continental erosionȦ l (index l for loss) Analogously, for the mass rate of change of the mantle water concentration, we havė Continental production depends both on the plate speed and on the rate of sediment subduction, which in turn depends on the continental area (see Section 4.1). A most simple parameterization of continental production that captures the key physics of the system is where the asterisk denotes parameters scaled to the present-day Earth, for examplė and where the index E denotes present-day Earth values. In Eq. A3, a fraction ξ 1 of the continental production rateȦ g is taken proportional to the continental area whereas the rest is taken proportional to the plate speed (derived in Section Appendix C). We account for the fact that subduction can only occur at those convergent plate boundaries that are not of continent-continent type by including the function L oo,oc (A) = L oo (A) + L oc (A), where L oo (A) and L oc (A) are the total lengths of ocean-ocean type subduction zones and ocean-continent type subduction zones, respectively (see Appendix B).
An essential part of continental erosion takes place on the surface. The rate of surface erosion increases with the surface area but depends, in addition, on the surface relief, the rate of rainfall, and the weatherability of the surface. Of these, we keep for simplicity the dependence on the surface area. There are other mechanisms of continental erosion of which we include subduction erosion for which the subducting plate scrapes off part of the overlying continent. The rate of subduction erosion mainly depends on the total length of ocean-continent type convergent plate boundaries. It would also depend on a friction coefficient but since we scale to the present Earth and assume that parameters are constant the friction coefficient will cancel. The subduction speed may matter but for plate speeds higher than the present one, the dependence is likely small (Stern, 2011). Thus, where ξ 2 is a constant.
We assume that the subducted water flux (regassing) to the mantle follows the same functional dependence on plate speed, sediment production and length of subduction zones as the continental production rate (see also Karlsen et al. 2019). This is a reasonable first-order approach, since the production rate of continental crust rock is taken directly proportional to the water subduction rate. The regassing rate should, of course, differ from the water subduction rate as water is consumed upon melting and lost through dehydration.
However, we assume for simplicity that the ratio between the two remains constant. The relative mass rates -scaled to their respective present-day Earth values -are, therefore, Mantle degassing, the release of water to the combined atmosphere-ocean reservoir depends, to first order, on the rate at which mantle material is transported into the melting region beneath mid-ocean ridges and other volcanic units and on the mantle water concentration.
We neglect degassing through volcanism in other geologic provinces.
where v * p is the plate speed scaled by its present-day value.
Equations 1 and 2 together with Eqs. A1 -A4 and A5 need values ofȦ l,E andȦ g,E andẇ l,E andẇ g,E for rescaling. Whereas the exact present-day rates may be of secondary importance, qualitative differences can result from either assuming an equilibrium or disequilibrium between present-day Earth sources and sinks. Even though a present-day equilibrium appears reasonable to assume, we cannot exclude that the sluggish Earth system only slowly converges towards an equilibrium state. We therefore introduce and note that a present-day Earth equilibrium implies f r = 1.
With Equations A1 -A8, the coupled continent-water system can qualitatively be described. A combination of these equations result in Eqs. 1 and 2. Particularly important are the constants ξ 1 and ξ 2 , which determine the strengths of positive and negative feedback, respectively.

Appendix B Length of subduction zones
Calculating the total length of subduction zones, we follow Höning et al. (2014) and assume that the total length of ocean-continent type subduction zones is proportional to the total length of continental margins (L oc ∼ L mar ). Modeling continents as randomly distributed spherical caps that may overlap, we find where L mar is the length of continental margins, R E is the radius of Earth. For n we use the number of continents on Earth of n = 7. Höning & Spohn (2016) have shown that L mar /L mar,max shown in Fig. B1 will not vary significantly with n for 5 ≤ n ≤ 9.
Eq. B1 has been derived as follows: Assume n spherical caps, each with area A 0 and circumference L 0 , Poisson randomly distributed on a sphere with surface area A surf . For a sufficiently large number n, the area covered by spherical caps A c is (Mecke & Stoyan, 2000;Schneider & Weil, 2008) where ν = n · A0 A surf , and the total length of margins L mar is Using the definition A = Ac A surf , Eq. B2 yields and combining Eq. B2 and Eq. B3 yields Furthermore, consider the radius of the base of the spherical cap r c and its height h c . The radius of the sphere R E is then Since L 0 = 2πr c and A 0 = 2πR E h c , we get (B7) Combining Eq. B5 and Eq. B7 yields and combining Eq. B4 and Eq. B8 results in Eq. B1.
With increasing continental area, continents increasingly overlap, which in turn reduces the total length of subduction zones. Following our assumption of L oc ∼ L mar , the total length of ocean-continent-type subduction zones has its maximum at a continental surface fraction of ≈ 0.37. For a smaller (larger) fraction we assume that the remaining length of subduction zones is of ocean-ocean (continent-continent) type (Fig. B1).
The derivative of Eq. B1 with respect to A is d dA

Appendix C Interior thermal evolution
The calculation of the thermal evolution of the mantle using Eq. 3 is based on the conventional thermal history models of terrestrial planets (e.g., Schubert et al., 2001) and includes terms describing the heat flow though the surface, heat production within the mantle, and the heat flow from the core.
For a constant (present-day Earth) continental crust coverage A E , the mantle heat production rate can be calculated following Korenaga (2008) where t is time in Gyr, Q 0 is the present-day Earth mantle heat production rate and q i and h i are the relative present-day heat production rates and half-lifes of the individual radiogenic elements 238 U, 235 U, 232 Th, and 40 K. Eq. C1 accounts for radiogenic decay but not for the transfer of radiogenic elements to the continental crust upon growth of the latter.
We account for the depletion of the mantle due to continental crust growth by writing where Q m (A, t) is the mantle radiogenic heat production rate as a function of time and continental coverage and where with Q C E the present-day heat production rate in the continental crust. In a scenario where we neglect the transfer of radiogenic elements, we set ξ 4 = 0.
The mantle heat flow is given by where T u and T s are the temperatures of the upper mantle and the surface, respectively, k is the thermal conductivity, A surf is the surface area of Earth, and δ u is the thickness of the upper boundary layer and given by where κ, g, α x , ρ, are the thermal diffusivity, gravitational acceleration, thermal expansivity, and density of the mantle, Ra c is the critical Rayleigh Number, ∆T is the temperature difference across the boundary layer, and η is the mantle viscosity given by where η 0 is a reference viscosity, E a and V a are the activation energy and volume, respectively, R g is the universal gas constant, T m and p m the temperature and pressure half-way through the mantle, T ref a reference temperature, and f w accounts for the dependence of the mantle viscosity on the concentration of water w. f w is given by where c 0 through c 4 are constants. Following boundary layer theory (e.g., Schubert et al., 2001), the plate speed is a function of the boundary layer thickness v p = 5.38κ d m δ 2 , where d m is the thickness of the mantle. Scaled to the present-day Earth, the mantle heat flow can, therefore, be written as where ∆T u = T u − T s and The parameterisation of the effects of thermal blanketing by the continental crust will complete the derivation of Eq. 3 and is introduced in Appendix D.

Appendix D Continental insulation
Considering thermal blanketing of the mantle by (partly) insulating continents, we begin by defining χ as the ratio of the present-day specific heat flow from the mantle through continental crust to that through oceanic crust where q o is the heat flow from the mantle through the oceanic crust and q c the heat flow from the mantle through the continental crust. We assume that χ remains constant throughout the evolution as it depends on the heat transfer properties of the oceanic and continental crust rather than on the thermal state of the mantle. We acknowledge, however, that χ may depend on the areal extent of the crust and the ratio between continental margins and old stable continental shield areas.
The heat flow from the mantle to the surface q s can then be written as The mantle heat flow q m,E has been defined by Eq. C10.
Since by definition or q o,E = ξ 3 q s,E and q c,E = (1 − ξ 3 )q s,E , we can write For the surface heat flow (Eq. D2), we get The present-day Earth surface heat flow is analogously: The relative heat surface heat flow is then (D7)

Appendix E Long-term carbon cycle
The plate speed determines the flux of mantle material into the melting region beneath mid-ocean ridges and thereby the mass rate of carbon outgassing from the mantle into the atmosphere where R man is the carbon reservoir in the mantle, V m the mantle volume, f deg the fraction of upwelling mantle material that degasses, L mor the total length of mid-ocean ridges, and d melt the depth below mid-ocean ridges where melting begins.
The weathering flux is calculated by (Foley, 2015) where P CO2 is the partial pressure of CO 2 in the atmosphere, P sat is the saturation pressure, E a,s is the activation energy of silicate weathering, T s is the surface temperature, f land is the land fraction, and F ws is the supply limit to weathering. Since we explicitly calculate surface erosion in this model by setting it proportional to the surface area of continental crust, we use f * land = A * . We note that this is a significant simplification in order to obtain first-order results without introducing additional equations. In reality, the emerged land fraction would also depend on additional parameters such as topography and surface water budget. The supply limit to weathering is given by where E max is the maximum physical erosion rate (in thickness per time), ψ is the fraction of reactable cations, m c is its average molar mass, and ρ r is the density of regolith (Foley, 2015).
Seafloor weathering is calculated as (Sleep & Zahnle, 2001;Foley, 2015) F Following classic long-term carbon cycle models (Sleep & Zahnle, 2001;Foley, 2015), the weathering fluxes transfer carbon from the atmosphere to the subduction zones considering half of that carbon is released during carbonate precipitation on the seafloor. Of the carbon that is subducted, we assume that another half is returned to the mantle while the rest is degassed at volcanic arcs back into the atmosphere.
The saturation pressure is calculated as where P sat0 and T s0 are reference values for the saturation pressure and surface temperature and m w and L w are the molar mass and latent heat of water, respectively. The surface temperature is calculated as T s = T s,E + 2(T e − T e,E ) + 4.6 P * CO2 0.346 − 1 , where T e is the effective temperature and given by where λ is the planetary albedo, σ is the Stefan-Boltzmann constant, and S is the solar irradiation, which is assumed to increase by a factor of 1/3 throughout 4.5 Gyr (Gough, 1981). Parameter values are given in Tables F1, F2, and F3.