A Nakayama

T Kodama

M Ikoma

Y Abe

Runaway climate cooling of ocean planets in the habitable zone: a consequence of seafloor weathering enhanced by melting of high-pressure ice

Keywords: planets and satellites, terrestrial planets -planets and satellites, oceans -planets and satellites, atmospheres c

Terrestrial planets covered globally with thick oceans (termed ocean planets) in the habitable zone were previously inferred to have extremely hot climates in most cases. This is because H 2 O high-pressure (HP) ice on the seafloor prevents chemical weathering and, thus, removal of atmospheric CO 2 . Previous studies, however, ignored melting of the HP ice and horizontal variation in heat flux from oceanic crusts. Here we examine whether high heat fluxes near the mid-ocean ridge melt the HP ice and thereby remove atmospheric CO 2 . We develop integrated climate models of an Earth-size ocean planet with plate tectonics for different ocean masses, which include the effects of HP ice melting, seafloor weathering, and the carbonate-silicate geochemical carbon cycle. We find that the heat flux near the mid-ocean ridge is high enough to melt the ice, enabling seafloor weathering. In contrast to the previous theoretical prediction, we show that climates of terrestrial planets with massive oceans lapse into extremely cold ones (or snowball states) with CO 2 -poor atmospheres. Such extremely cold climates are achieved mainly because the HP ice melting fixes seafloor temperature at the melting temperature, thereby keeping a high weathering flux regardless of surface temperature. We estimate that ocean planets with oceans several tens of the Earth's ocean mass no longer maintain temperate climates. These results suggest that terrestrial planets with extremely cold climates exist even in the habitable zone beyond the solar system, given the frequency of water-rich planets predicted by planet formation theories.

INTRODUCTION

The Earth's climate system is generally thought to be stabilized by a carbonate-silicate geochemical cycle of carbon (hereafter called the carbon cycle). On geological timescales, the amount of the greenhouse gas CO 2 is determined by a balance between degassing flux through volcanism and sinking flux through chemical weathering. Since chemical weathering becomes more efficient with temperature, a negativefeedback mechanism operates to keep the CO 2 partial pressure at low levels and, thus, to maintain the temperate cli-⋆ E-mail: anakayama@eps.s.u-tokyo.ac.jp (AN) mate (Walker et al. 1981). In the present Earth, weathering occurs mainly on continents [START_REF] Caldeira | [END_REF]. Beyond the solar system, however, there must be continent-free terrestrial planets completely covered with oceans in the habitable zone. In this study we refer to the conventional habitable zone defined based on 1-D radiativeconvective models, ranging from 0.34 to 1.06 times the present solar insolation at the Earth's orbit (Kasting et al. 1993;Kopparapu et al. 2013), as the habitable zone. Given diverse water-supply processes and their stochastic nature, terrestrial exoplanets must be diverse in ocean mass. Indeed, many recent theories of planet formation predict that terrestrial exoplanets could have much more water than the Earth (see recent reviews by O'Brien et al. ( 2018) and Ikoma et al. (2018)). N-body simulations of late-stage terrestrial planet accretion including the supply of water-rich planetesimals beyond the snowline demonstrate that terrestrial planets with oceans of ten to several hundred Earth's ocean masses might be common in the habitable zone (Raymond et al. 2004(Raymond et al. , 2007)). On the Earth, there would be no lands if the ocean mass were three times larger than the present (i.e., 0.023 % of the Earth's mass) (Maruyama et al. 2013;Kodama et al. 2018).

What is the climate like on terrestrial planets completely covered with oceans and what influence does ocean mass have on the climate? Such terrestrial planets are called ocean planets, hereafter, whereas ones covered partially with oceans like the Earth are called partial ocean planets (Kuchner 2003;Léger et al. 2004). On ocean planets, seafloor weathering, instead of continental weathering, would control the planetary climate. The role of seafloor weathering in climate is of interest in this study. In particular, we focus on the influence of high-pressure (HP) ice of H 2 O such as ice VI and VII on the seafloor weathering.

Planets with larger water amounts than a certain threshold have the HP ice on the seafloor, provided the ocean has a steady, isothermal or adiabatic structure (see Fig. 2 of Léger et al. 2004). Since the HP ice is a solid heavier than its counterpart liquid, a solid layer is formed between the ocean and oceanic crust and prevents seafloor weathering (Alibert 2014;Kaltenegger et al. 2013).

Climates of ocean planets without geochemical interaction and carbon cycle between the ocean-atmosphere system and silicate mantle were previously investigated. Wordsworth & Pierrehumbert (2013) and Kitzmann et al. (2015) explored the effect of dissolution of CO 2 into a cationpoor ocean and found that the CO 2 pressure decreases with increasing temperature for a given carbon inventory in the atmosphere-ocean system (see Fig. 2 of Kitzmann et al. (2015)). This suggests that such a climate system is an unstable one with a positive feedback cycle. Kite & Ford (2018) considered supply of cations to the ocean, which strongly affects ocean chemistry, in the initial, hot stage after solidification of the magma ocean. They showed that large cation concentration enhances the positive feedback and leads to destabilizing planetary climate into hot one for a large CO 2 inventory (∼ 100 bars) in the atmosphere-ocean system even for stellar insolation comparable to the present Earth.

However, whether the layer of HP ice really exists and prevents seafloor weathering completely must be verified through a detailed consideration of heat transfer and rheology in the HP ice layer. Noack et al. (2016) examined the stability of the HP ice layer by performing non-steady, onedimensional simulations of heat transfer, including the melting of HP ice, in the layer (liquid H 2 O + HP ice) above the oceanic crust (collectively called the H 2 O layer, hereafter). They found that the heat flux from the oceanic crust is too high for steady heat transport in the HP ice and, thus, the heat is temporarily stored near the bottom of the H 2 O layer, which results in melting the HP ice. Since the resultant melt is lighter than its surroundings, an upwelling flow of partially molten HP ice occurs. Such a possibility has been investigated also in studies of large icy moons in the solar system, in particular, Ganymede, which propose that solid and liquid coexist via melt production within the HP ice layer, bringing about a melt-buoyancy-driven upwelling flow in the interior.

To evaluate the efficiency of heat transport by the meltbuoyancy-driven flow, Choblet et al. (2017) performed 3-D simulations of thermal convection in the HP ice layer, including the effect of melting of the HP ice. In their simulations, they assumed and mimicked a permeable flow in the HP ice by extracting the generated melt instantaneously to the above ocean. Then, they demonstrated that melt is mostly generated on the oceanic crust and the permeable flow dominates the heat transport. Recently, Kalousová et al. (2018) performed 2-D convection simulation of a water-ice mixture to investigate the behavior of the generated melt in the HP ice layer. They demonstrated that heat is efficiently transported by the melt-buoyancy-driven convective and permeable flows and water is exchanged throughout the HP ice layer. In this study, we call those flows the sorbet flow, since they are flows of a water-ice mixture. The sorbet flow occurs for the small thickness of the HP ice ( 200 km) and large heat flow ( 20 mW m -2 ) for Ganymede-like icy bodies. Nusselt-Rayleigh number scaling supports that such a sorbet flow likely occurs also for ocean planets with Earthlike geothermal heats (80 mW m -2 in the present Earth's mean mantle heat flow) and thicker HP ice. Hence, seafloor weathering likely occurs for ocean planets with the HP ice.

Horizontal variation is another important effect ignored previously. In particular, for planets where plate tectonics works, the heat flow from oceanic crusts is highest at mid-ocean ridges and decreases with distance from there. The heat flow near mid-ocean ridges can be high enough to melt the HP ice. Then, the seafloor temperature is fixed close to the melting temperature for the pressure at the seafloor (hereafter, the seafloor pressure). This temperature is much higher than one obtained from inward integration of the adiabat from the oceanic surface to the seafloor. Higher seafloor temperature results in more efficient seafloor weathering, according to the temperature dependence of seafloor weathering inferred based on dissolution experiments of basalt (Brady & Gíslason 1997;Gudbrandsson et al. 2011) and geological evidence (Coogan & Dosso 2015;Krissansen-Totton & Catling 2017). Hence, the seafloor weathering can remove atmospheric CO 2 efficiently, provided such a molten region is sufficiently wide.

This study is aimed at evaluating the role of the HP ice in seafloor weathering and climate for ocean planets with a focus on the effects of the liquid-solid coexistence region maintained by the sorbet flow and the horizontal variation in heat flux from the oceanic crust. The rest of this paper is organized as follows: In section 2, we describe our model to simulate the ocean layer structure and planetary climate. In section 3, we show the behavior of the HP ice with a focus on the area where melting occurs. In section 4, we show the impacts of seafloor weathering with the HP ice on the planetary climate. In section 5, we discuss surface environments of ocean planets, caveats of the model, and implication of our results for terrestrial exoplanets. In section 6, we conclude this study.

CLIMATE MODEL

We consider an Earth-size ocean planet with various amounts of H 2 O and CO 2 . Of special interest in this study is the impact of ocean mass, M oc , on the planetary climate including surface temperature, T s , and CO 2 partial pressure, P CO 2 . We assume that the planet is almost Earth-like, namely, a terrestrial planet with the Earth's mass and internal composition orbiting at 1 AU far from a Sun-like star, except for the ocean mass. Our climate model consists of four components: (1) internal structure integration that determines the thickness of the HP ice layer (section 2.1);

(2) seafloor environment modeling that determines the area where seafloor weathering works when the HP ice is present (section 2.2); (3) carbon cycle modeling that calculates P CO 2 (section 2.3); (4) atmospheric modeling that calculates T s (section 2.4).

Ocean structure model

The hydrostatic structure of the ocean is determined by

dP dr = -gρ, (1)
dm dr = 4πr 2 ρ, ( 2
)
where r is the radial distance from the planetary center, P and ρ are the pressure and density, respectively, m is the cumulative mass, and g is the gravity (g = Gm/r 2 ; G being the gravitational constant). Its thermal structure is assumed to be adiabatic:

dT dr = - αgT C P , (3)
where T is the temperature, α is the thermal expansivity, and C P is the heat capacity. Two of the three boundary conditions are T = T s and P = P s at r = R p , where P s is the surface pressure and R p is the planet radius. Here P s is the sum of the background pressure (P n = 1 × 10 5 Pa) and vapor pressure P H 2 O for T s , which is taken from Nakajima et al. (1992); P CO 2 is relatively small. The inner boundary condition is m = 0 at r = 0. This means that we must continue the integration until the planet's center, although we are interested only in the ocean layer. The planet consists of a H 2 O ocean layer, a rocky mantle, and an iron core. Regarding the equations of state for the materials and chemical phases, we mostly follow Valencia et al. (2007a). The details are given in Appendix A.

The red line of Fig. 1 shows the calculated relationship between the surface temperature and critical ocean mass beyond which HP ice exists (see section 2.5 for the numerical procedure), which is abbreviated to COM-HP hereafter. It turns out that HP ice exists for an Earth-like planet with an ocean of more than ∼20 to ∼100 M oc,⊕ , depending on surface temperature. To see the sensitivity to the thermal structure of the ocean, we also show the result for an isothermal ocean. The difference in COM-HP between the isothermal and adiabatic cases is ∼ 1-30M oc,⊕ for T s = 280-400 K. Even for the two extreme cases, the difference is small enough not to change our conclusions.

Seafloor environment model

Near a mid-ocean ridge, heat flow from below is so high that the HP ice would be incapable of transporting the heat by thermal conduction nor convection and consequently become molten. If liquid water exists together with ice, the Surface temperature [K] Figure 1. The critical ocean mass (COM-HP) in the unit of the Earth's ocean mass (M oc,⊕ ), beyond which high-pressure (HP) ice appears deep in the ocean is shown as a function of surface temperature (red line). Note that the result for an isothermal liquid ocean is also shown (blue line) to confirm that our calculation reproduces the result of Kitzmann et al. (2015) well.

heat can be transported efficiently by a sorbet flow, as described in Introduction. The HP ice far from a mid-ocean ridge remains solid because of low heat flow. Hence, there is a critical distance beyond which or a critical heat flow, q cr , below which the HP ice remains solid.

A schematic illustration of our seafloor environment model is shown in Fig. 2. Here we assume that (1) the heat transport is steady and vertically one dimensional, (2) the composition and phase of H 2 O are vertically homogeneous in the HP ice region, and (3) the sorbet flow dominates the heat transport in the solid-liquid coexistence region (called the sorbet region, hereafter), whereas solid-state convection or conduction occurs in the HP ice region. These assumptions are consistent with results of previous hydrodynamical simulations (Choblet et al. 2017;Kalousová et al. 2018). We discuss their validity and impact on our conclusion in section 5.3.1.

Critical heat flow

First, we determine the critical distance or critical heat flow. Namely, according to its definition, we find the point at which convection nor conduction can hardly transport the heat inside the HP ice region. Figure 2b shows qualitative temperature profiles in the HP ice: At the critical distance, since the ice-liquid mixture on the oceanic crust is in phase equilibrium, the temperature is equal to the melting temperature. Also, the temperature at the top of the HP layer is the melting one, by definition.

To determine the thermal structure of the HP ice region, we adopt a similar approach with that used by Fu et al. (2010) who investigated the structure of the icy mantle of an ocean planet with a frozen surface, although they ignored horizontal variation in heat flux. Unlike Fu et al. (2010), we take into account the case where the HP ice layer is wholly conductive, ignore the upper thermal boundary layer, and in the infinitesimally thin layer on the side of the high-pressure ice at the boundary between the "sorbet" and "high-pressure ice" regions.

In panel (a), arrows represent the direction and dominant mechanisms of heat transport. In the sorbet region, ice and liquid coexist and thus the temperature is fixed at the melting point of H 2 O. In panel (b), Ra and Ra cr represent the Rayleigh number and the critical Rayleigh number, respectively. The red and orange solid lines represent thermal structure when the high-pressure ice layer is convective and conductive, respectively. The blue and black solid lines represent the melting line of H 2 O and adiabatic thermal structure of the liquid ocean, respectively.

consider the different boundary condition for the bottom of the HP layer. The details are described below.

The mechanism of heat transport depends on the Rayleigh number, Ra, which is defined as [START_REF] Turcotte | Geodynamics -2nd Edition[END_REF])

Ra = gαρD 3 ∆T HP κη , ( 4
)
where D is the thickness of the HP ice layer, κ is the coefficient of thermal diffusivity, η is the viscosity, ∆T HP = T mel BB -T mel TB and T mel TB and T mel BB are the melting point temperatures for the pressures at the top and bottom of the HP ice layer, respectively. For the melting point temperature T mel , we use the formula from Dunaeva et al. (2010),

T mel = a 1 + a 2 P + a 3 ln P + a 4 P -1 + a 5 √ P, (5)
where P is the pressure in bar and the values of coefficients are summarized in Table 1. We assume that the phase transition from ice VI to VII occurs at the triple point of liquid/ice VI/ice VII, the pressure of which is 22160 bars. The thermal diffusivity is defined by κ = k/ρC P , where k is the thermal conductivity. For C P of ice VI and VII, we use the expression derived by Fei et al. (1993). For k, we adopt a constant value of 3.8 Wm -1 K -1 , which is its typical value for ice VII under 2.5 GPa and 300 K (Chen et al. 2011), for simplicity. For η of ice VII, which is poorly constrained, we adopt a dislocation model for the viscosity of phase VI, which is the highest phase of the HP ice measured so far (Durham et al. 1997) :

η(P η , T η ) = Bζ -3.5 exp (E * + P η V * ) RT η , (6)
where B (= 6.7 × 10 19 Pa 4.5 s) is a constant, ζ (= 2.0 × 10 6 Pa) is a characteristic shear stress (Fu et al. 2010), R is the ideal gas constant, E * (= 110 kJ mol -1 ) and V * (= 1.1 × 10 -5 m 3 mol -1 ) are the activation energy and volume (Durham et al. 1997), respectively, and T η and P η are the temperature and pressure at deformation, respectively.

Because the viscosity contrast in the HP ice layer is relatively small, the small viscosity contrast prescription can be used (Fu et al. 2010). For P η and, T η , we use the averaged values for the HP ice layer (Dumoulin et al. 1999). In this study, we assume the value of the critical Rayleigh number, Ra cr , is 2000. When Ra < Ra cr , conduction dominates heat transport and, thus, q cr is given as

q cr = k ∆T HP D . (7)
When Ra > Ra cr , since convection occurs, we assume the adiabatic temperature gradient (i.e., Eq.[3]). Near physical boundaries, however, since convective motion is prevented, thermal boundary layers are formed, where conduction transports heat. In this study, we consider the presence of a boundary layer only on the bottom of the HP ice layer (BBL), where the temperature gradient is given by

dT dr = - q k (8)
and q is the heat flux. Integrating Eq. (3) inwards from the top of the HP ice layer and Eq. ( 8) outwards from the surface of the oceanic crust, we determine the BBL's thickness, δ , and the temperature difference in the BBL, ∆T BBL at the crossover point for a given q (see Fig. 2b).

Given that the BBL is marginally stable against convection, Ra = Ra cr in the BBL, namely, Ra cr = gαρδ 3 ∆T BBL /κη BBL , which comes to be

δ = κη BBL Ra cr gαρ∆T BBL 1/3 , (9)
where η BBL is the viscosity of the HP ice in the BBL and calculated with the intermediate values of temperature and pressure between the top and bottom of the BBL. If the set of δ and ∆T BBL for a given q satisfies Eq. ( 9), the value of q corresponds to q cr , which is also written as

q cr = k ∆T BBL δ = kκRa cr gαρ • η BBL δ 4 . ( 10
)
Note that Fu et al. (2010) considered a boundary layer under the top of the HP ice layer in addition to BBL. We discuss the difference in temperature structure in the HP ice layer between this study and Fu et al. (2010) and its impacts on our conclusion in section 5.3.1.

Effective weathering area

Same as in the Earth, the oceanic crust is assumed to form via eruption of hot mantle rock only at the mid-ocean ridge.

As it moves away from the mid-ocean ridge toward the trench, the oceanic crust is cooled by seawater. Here we define a non-dimensional effective weathering area, f oc , as the area of the sorbet region (i.e., q > q cr ) relative to the whole area of oceanic crust. A constant rate of oceanic crust production being assumed, f oc is equivalent to the ratio of the period during which q ≥ q cr to the residence time of the oceanic crust, τ.

To calculate f oc , we model the cooling of the oceanic crust, adopting the semi-infinite half-space cooling model [START_REF] Turcotte | Geodynamics -2nd Edition[END_REF]: This model assumes that the crust cools only by vertical heat conduction. The heat flux from the oceanic crust is given by [START_REF] Turcotte | Geodynamics -2nd Edition[END_REF]

q(t) = k rock (T sol -T floor ) √ πκ rock t ≡ A √ t , ( 11
)
where t is time, k rock (= 3.3 W m -1 K -1 ) and κ rock (= 1.0 × 10 -6 m 2 s -1 ) are the thermal conductivity and thermal diffusivity of the oceanic crust, respectively, T floor is the seafloor temperature, and T sol is the potential temperature of the mantle, for which we assume the peridotite dry solidus at the seafloor pressure, which was parameterized by Hirschmann et al. (2009). This assumption is made just for simplicity. The influence of the assumption on planetary climate is discussed in section 5.2. From Eq. (11), the length of time required for q to decrease to q cr , which is denoted by t cr , is given by t cr = A 2 /q 2 cr . Also, if the mean mantle heat flow, q, is defined as q ≡ τ -1 τ 0 q dt, the residence time τ is given as a function of q as τ = 4A 2 / q2 . Thus, the effective weathering area is given as

f oc ≡ t cr τ = 1 4 q q cr 2 . ( 12
)
In some cases, calculated t cr happens to be larger than τ, which means the oceanic crust is fully covered with the solidliquid mixture (i.e., the sorbet). In such cases, we set f oc = 1.

From Eq. ( 12), it turns out that when q cr > q/2, solid HP ice appears near the trench. In the next section, we use f oc in the carbon cycle model. In this study, the mean mantle heat flow q is a free parameter. As the fiducial value, we use q = 80 mW m -2 , which is the value for the present Earth. Note that q cr is independent of q, according to Eq. ( 12).

Carbon cycle model

In order to investigate planetary climate, we develop a carbon cycle model by modifying the Earth's carbon cycle model of [START_REF] Tajika | [END_REF]. Since we focus on continent-free terrestrial planets, we add the effect of seafloor weathering and neglect the continental reservoir of carbon and the effect of continental weathering. In addition, we consider the presence of the HP ice and pressuredependent degassing. Same as [START_REF] Tajika | [END_REF] and Sleep & Zahnle (2001), we perform box-model calculations of carbon circulation among reservoirs and find the equilibrium states.

Carbon reservoirs

We consider four reservoirs, which include the atmosphere, ocean (liquid water plus HP ice), oceanic crust (basalt), and mantle. Between the atmosphere and ocean, however, the carbon partition is assumed to be always in equilibrium, which is described in detail in Appendix B. The equilibrium value of the CO 2 partial pressure P CO 2 depends on the number of cations dissolved in the ocean (e.g., [START_REF] Zeebe | CO2 in seawater: equilibrium, kinetics, isotopes[END_REF], for which we assume the present Earth's value, although the supply of cations via continental weathering never occurs in ocean planets. We have confirmed that overall results are insensitive to the number of dissolved cations (even in the case with no cations in the ocean). This is because the ocean reservoir is much small relative to the whole planetary carbon reservoir.

Carbon dissolved in the ocean is carried to the seafloor in the form of CO 2 ice, into which aqueous CO 2 is converted in the sorbet region (Bollengier et al. 2013). We assume that the carbon circulation in the sorbet region occurs quickly enough that it never affects the mass balance and also planetary climate. Detailed discussion of the CO 2 circulation is given in section 5.3.1.

The origin of volatiles in terrestrial planets has been highly debated so far, even for the Earth (e.g., O'Brien et al. 2018). Since possible candidates such as carbonaceous chondrites and comets include both carbon compounds and water, we assume that the total mole number of carbon contained in the whole planet, C total , is proportional to ocean mass, namely

C total = γ n oc,⊕ M oc M oc,⊕ , ( 13
)
where γ is the CO 2 /H 2 O molar ratio in the source of volatiles of the planet, n oc,⊕ (= 7.6 × 10 22 mol) is the molar quantity of H 2 O in the Earth ocean mass, and M oc,⊕ (= 1.37 × 10 21 kg) is the Earth's ocean mass. Using the data and estimation published, we can estimate γ is to be 0.22 for carbonaceous chondrites (Jarosewich 1990), 0.71 for comae of comets (Marty et al. 2016), and 0.19 for Earth composition [START_REF] Tajika | [END_REF]. We use the Earth-like value (γ = 0.19) as the fiducial value. Dependence of planetary climate on γ is discussed in sections 4.

Carbon budget

The mass balance among those reservoirs is expressed as

d(C atm +C oc ) dt = F D + F M -F SW , (14)
dC bs dt = F SW -F R -F M , ( 15
)
dC man dt = F R -F D , (16)
C total = C atm +C oc +C bs +C man , (17)
where C atm , C oc , C bs , and C man are the mole numbers of carbon contained in the atmosphere, ocean, oceanic basalt, and mantle, respectively, and F SW , F D , F R , F M are the carbon fluxes due to seafloor weathering, degassing from the mid-ocean ridge, regassing via subduction into mantle and metamorphism that leads to degassing from volcanic arc, respectively. Those equations are solved for a given value of C total . We adopt the degassing, regassing, and metamorphism models from [START_REF] Tajika | [END_REF], where each flux is expressed as

F D = K D A S C man , (18)
F R = β τ C bs , (19)
F M = 1 -β τ C bs . ( 20
)
Here K D is the molar fraction of carbon degassing as CO 2 from the erupting magma per unit area. We take into account the dependence of K D on seafloor pressure (i.e., ocean mass), the detail of which is described in Appendix C. β is the regassing ratio defined as the molar fraction of carbonate regassed into the mantle in the total subducting carbonate. We adopt the present Earth's value of β (= 0.4) estimated by [START_REF] Tajika | [END_REF]. A S is the seafloor spreading rate, which is simply given by

A S = A 0 τ , (21)
where A 0 is the whole area of the seafloor. We assume that A 0 is the present Earth's value (= 3.1 × 10 14 m 2 ) from McGovern & Schubert (1989) and calculate τ from the relation τ = 4A 2 / q2 for a given q.

The seafloor weathering rate F SW depends on seafloor temperature T floor as (Brady & Gíslason 1997)

F SW = F * SW f oc exp E a R 1 T 0 - 1 T floor , (22)
where F * SW is the present Earth's seafloor weathering rate, f oc is the effective weathering area given by Eq. ( 12), E a is the activation energy, and T 0 (= 289 K) is the reference seafloor temperature that corresponds to the surface temperature obtained by the atmospheric model with the present Earth's condition. F * SW estimated from deep-sea cores is 1.5-2.9 × 10 12 mol yr -1 (Alt & Teagle 1999; [START_REF] Staudigel | [END_REF]Gillis & Coogan 2011). In this study, we use F * SW = 2.0 × 10 12 mol yr -1 .

The activation energy E a is uncertain and its reported value ranges between 30 and 92 kJ mol -1 . Brady & Gíslason (1997) firstly determined E a experimentally to be 41 kJ mol -1 . Recent inversion methods using geological evidence support a relatively high value of E a : Precisely, strontium and oxygen isotopes in carbonates indicated E a = 92 ± 7 kJ mol -1 (Coogan & Dosso 2015). Also, several proxies reflecting the surface and seafloor temperatures, atmospheric CO 2 , and oceanic pH showed E a = 75 +22

-21 kJ mol -1 (Krissansen-Totton & Catling 2017). Those values are also consistent with estimates from laboratory experiments for the dominant minerals in the oceanic crust [START_REF] Brantley | Reaction Kinetics of Primary Rock-Forming Minerals under Ambient Conditions[END_REF]. In contrast, an experimental study of basalt dissolution in the moderate pH range reported the relatively small E a of 30 kJ mol -1 (Gudbrandsson et al. 2011). In this study, we use E a = 41 kJ mol -1 as the fiducial value according to previous studies (e.g., Foley 2015) and vary it over the range between 30 and 92 kJ mol -1 . We ignore the pH dependence of seafloor weathering since it is known to be small in the pH range between 4 and 10 (Gudbrandsson et al. 2011).

The seafloor temperature also depends on the surface temperature, T s , because we assume that the temperature structure of the ocean is adiabatic (see also § 2.1). We calculate T s as a function of P CO 2 , as described in detail in section 2.4. On the area of the seafloor beneath the sorbet region, T floor is equal to the melting temperature at the seafloor pressure.

Atmospheric model

In this study, we use the open-source code for 1-D radiativeconvective climate models, Atmos 1 , developed by Kasting and his collaborators (Kasting et al. 1993;Kopparapu et al. 2013;Ramirez et al. 2014). This code calculates radiative fluxes in vertically spacing layers of the atmosphere, using the two-stream approximation with the coefficients for radiative absorption and scattering by gaseous molecules updated by Kopparapu et al. (2013). We assume a 1-bar N 2 atmosphere with various partial pressures of CO 2 . The distribution of the relative humidity of water vapor is treated according to the empirical Manabe-Wetherald model which assumes the surface relative humidity of 0.8, based on the present Earth's atmosphere (Manabe & Wetherald 1967;Pavlov et al. 2000). According to Kopparapu et al. (2013), we use the surface albedo of 0.32, which implicitly includes the effects of present-day Earth water clouds. We use the present insolation flux at the Earth's orbit S ⊙ (=1360 W m -2 ) and the present Sun's spectrum as the fiducial value and spectrum model, respectively. The other model settings are the same as those adopted in Ramirez et al. (2014). Then, we calculate equilibrium values of T s as a function of P CO 2 for given stellar insolation, using a time-stepping approach with moist convective adjustment (Pavlov et al. 2000). We have confirmed that our calculated T s is almost the same with sufficient accuracy as that from Ramirez et al. (2014). We discuss the uncertainties and impacts of stellar insolation, surface albedo, and relative humidity in sections 4.2, 5.3.2, and 5.4.

Numerical procedure

In summary, for given values of ocean mass M oc and mean mantle heat flow q, we determine the climate of the ocean planet by the following procedure.

(i) For trial values of surface temperature T s and surface pressure P s , we integrate Eqs. ( 1)-( 3) inward from the surface to determine temperature as a function of pressure in the ocean (see § 2.1). We find a level where the adiabat crosses the melting temperature of ice. The layer between the crossover level and the oceanic crust surface consists of HP ice. Then, the seafloor pressure P floor and the thickness of the HP ice layer D are determined. If the adiabat reaches the oceanic crust surface before crossing the ice melting curve, the planet has no ice in the deep ocean. The numerical integration is performed with a 4th-order Runge-Kutta method. The size of the interval is chosen so that the pressure at the crossover point is determined with < 0.1 % accuracy.

(ii) When the HP ice is present, from the seafloor environment model, we determines the critical heat flow q cr (or the area of the sorbet region) from Eq. ( 7) or ( 10), depending on Ra ( § 2.2). Then, we obtain the effective weathering area f oc by substituting q cr and q in Eq. ( 12). Also, we obtain the seafloor temperature T floor in the sorbet region by substituting P floor in Eq. ( 5).

(iii) In the carbon cycle model ( § 2.3), using T floor and f oc obtained above, we perform a time integration of Eqs. ( 14)-( 16) and determine the carbon partition among the atmosphere, ocean, oceanic crust, and mantle. Then, from the calculated P CO 2 , we obtain a new value of T s (and thereby P s ) from the atmospheric model ( § 2.4). If the new value of T s differs by > 0.01 K from the trial value of T s , we return to Step (i) and repeat the above procedures with the new T s . The time integration is performed with a Euler method and the interval size is chosen so that the time difference in the molar number of carbons is smaller than 0.1 % for all the reservoirs.

(iv) Once all the time derivatives in Eqs. ( 14)-( 16) become zero, we judge the solution as an equilibrium state. If the surface temperature drops below 273 K, we also stop the time integration and regard the solution as a snowball state.

We start time-stepping calculations at arbitrarily high P CO 2 (i.e., in a warm condition) for finding equilibrium solutions. We have confirmed that the results are insensitive to choice of the initial condition, provided a sufficiently high CO 2 pressure (P CO 2 > 10 bars) is adopted. (The carbon cycle and climate stability in the snowball state are discussed in section 5.3.3.) In most of our simulations, response against perturbations for the carbon budget in the atmosphereocean system is mainly controlled by regassing, the timescale of which is ∼ τ/β = 250 Myr for q = 80 mW m -2 . Thus, an equilibrium state is achieved on a timescale of the order of Gyr, which is also consistent with results shown in Foley (2015).

The parameters and constants with their values adopted in this study are summarized in Tables 2 and3, respectively. The upper limit for ocean mass M oc that we consider is 200 M oc,⊕ . The reasoning is as follows: We suppose that plate tectonics is working on the planet. Although still not fully understood, an increase in water has negative effects on plate tectonics. In particular, it leads to reducing crustal production and degassing, since the solidus temperature of the mantle material increases with pressure (Kite et al. 2009;Noack et al. 2016). According to Noack et al. (2016), crustal production completely ceases for an Earth-mass planet with the ocean layer thicker than approximately 400 km, if plate tectonics operates. The ocean mass of 200 M oc,⊕ that we adopt here corresponds to the ocean layer of ∼350 km for T s = 300 K. We do not consider ocean planets with more massive oceans because such planets are expected to have no geochemical cycle. For planetary climates with no geochemical cycle, see Kitzmann et al. (2015) and Kite & Ford (2018).

Note that we assume a spherically symmetric structure in the internal structure modeling, while we consider the presence of the sorbet and HP ice regions in the deep ocean in the seafloor environment modeling. Such selfcontradiction, however, has little influence on our whole modeling. This is because only the thermal structure above the HP ice layer is of interest in this study and the equations of state of water, rock, and iron are rather insensitive to temperature.

MELTING OF THE HP ICE

We first investigate the behavior of the HP ice with a focus on the effective weathering area, which is a controlling factor for seafloor weathering. Here we do not use the carbon cycle model, but, instead, perform calculations for fixed values of the surface temperature T s . Figure 3 shows the calculated thickness of the HP ice layer D (left column), the critical heat flow q cr (middle column) and effective weathering area f oc (right column) as a function of ocean mass M oc for T s = 300 K (top) and as a function of T s for M oc = 200M oc,⊕ (bottom). In those calculations, the mean mantle heat flow q is assumed to be 80 mW m -2 .

Dependence on Ocean Mass

The overall dependence on ocean mass is as follows. As shown in Fig. 3a, the HP ice is present, if M oc 45M oc,⊕ . Its thickness increases almost linearly with ocean mass and reaches 247 km at M oc = 200M oc,⊕ . In Fig. 3b, the critical heat flow is found to be zero for M oc 45M oc,⊕ , because of no HP ice, and then increase with ocean mass, up to about 80 mW m -2 (≃ q) at M oc = 200M oc,⊕ . In Fig. 3c, the effective weathering area is found to be unity until M oc ≃ 139M oc,⊕ and rapidly decrease to about 0.2 at M oc = 200M oc,⊕ .

A jump in q cr is found at M oc ≃ 74M oc,⊕ in Fig. 3b. At that point, the heat transport mechanism in the HP ice above the critical point (i.e., q = q cr ) changes from conduction to convection. For M oc 74M oc,⊕ (D ≃ 55 km), the HP ice layer is thin enough and, therefore, the temperature difference ∆T HP (= T mel BB -T mel TB ) is small enough for conduction to transport the heat flux from the oceanic crust. However, as shown in Fig 3c, meaning that the HP ice is entirely molten (i.e., f oc = 1), that is, the seafloor is covered entirely with the sorbet for M oc 74M oc,⊕ (see the text just below Eq. [12]). Note that discontinuities in q cr or dq cr /dM oc found at M oc ≃ 86 and 103M oc,⊕ come from those in the melting curve of H 2 O at the phase boundaries of ice VI/VII. The critical heat flow exceeds q/2 at M oc ≃ 139 M oc,⊕ (D ≃ 160 km), until which the effective weathering area is unity, and then increases further with ocean mass. Such an increase in q cr occurs because the Rayleigh number in the HP ice layer increases. Thus, the effective weathering area decreases with ocean mass, but never becomes zero until M oc = 200M oc,⊕ . This means that water-rock reactions between water and rock including the seafloor weathering are possible, despite the presence of the thick HP ice, because the sorbet region also exists near the mid-ocean ridge.

q cr ≃ 10 mW m -2 < q/2 at M oc 74M oc,⊕ ,

Dependence on Surface Temperature

The three lower panels of Fig. 3 show the dependence on the surface temperature for M oc = 200M oc,⊕ . The HP ice thickness decreases, as the surface temperature increases, as shown in Fig. 3d. At T s ≃ 320 K, the curve is a bit inflected. This is due to the phase change of HP ice from ice VI to ice VII.

In Figs. 3e and 3 f , we find a maximum of the critical heat flow and a minimum of the effective weathering area, respectively, at T s ≃ 390 K. As indicated in Eq. ( 10), q cr depends on η BBL and δ , both of which decrease with T s . For T s 390 K, δ 4 decreases more rapidly than η BBL and, thus, q cr increases with T s . In contrast, for T s 390 K, the latter dominates over the former, so that q cr decreases. At T s ≃ 390 K, ∂ (η BBL /δ 4 )/∂ T s = 0. The behavior of the effective weathering area can be readily understood from Eq. ( 12), namely, f oc ∝ q -2 cr . The minimum is f oc ≃ 0.04. In Fig. 4, we show the maximum of critical heat flow q cr,max and minimum of effective weathering area f oc,min as a function of ocean mass for q = 80 mW m -2 . Here we show only the results for the case of convective HP ice for M oc > 100M oc,⊕ because the critical heat flow due to conduction is small. While q cr,max is found to monotonically increase with M oc , f oc,min begins to drop with M oc at M oc ≃ 128M oc,⊕ , which is smaller than in the case of T s = 300 K because of difference in T s . The blue line in Fig. 4 indicates that even the minimum of f oc is unity for M oc 128M oc,⊕ , which means that the HP ice is entirely molten and the seafloor is completely covered with the sorbet, regardless of surface temperature, in such an ocean mass range for the Earth-like mean mantle heat flow ( q = 80 mW m -2 ). Also, f oc,min > 0, meaning that seafloor weathering works, even if M oc = 200M oc,⊕ .

SEAFLOOR WEATHERING ENHANCED BY THE HP ICE

Consequence of carbon cycle

Here we examine the planetary climate based on the carbon cycle including the effective weathering area obtained above.

The calculation results for q = 80 mW m -2 , γ = 0.19, and S = S ⊙ and 0.9S ⊙ are shown in Fig. 5, where (a) the surface and seafloor temperatures, (b) the seafloor weathering flux and effective weathering area, and (c) the partial pressure of atmospheric CO 2 are plotted as functions of the ocean mass. In Fig. 5a, two obviously different states are found: One is the state with T s > 273 K, where the carbon cycle is in a steady state, the other, as indicated by a shaded area, is the state with T s = 273 K, where the carbon cycle calculation is artificially stopped at T s = 273 K because the surface ice is expected to form (see also § 2.5). The former is called the equilibrium state and the latter is called the snowball state in this study. In this case, the HP ice begins to form at M oc = 86M oc,⊕ . It turns out that the formation of the HP ice has a drastic effect on the carbon cycle and determines which state is achieved.

In the case of no HP ice (i.e., M oc < 86M oc,⊕ ), both the surface temperature and CO 2 partial pressure increase with ocean mass. An equilibrium state is achieved for a given ocean mass via a negative feedback loop such that an increase in P CO 2 raises the surface temperature, which leads to a rise in seafloor temperature, which enhances seafloor weathering flux, which finally reduces the atmospheric CO 2 . The larger the ocean mass, the larger the total carbon inventory C total is (see Eq. [13]). Since an increase in C total enhances the degassing flux of CO 2 (see Eq. [18]), the surface temperature consequently raises with ocean mass. This is, in other words, because the enhancement of the degassing flux dominates over the increase in seafloor temperature in the case of γ = 0.19. While the outcome depends on γ, we have confirmed that this trend is the same also in the case of one-tenth of the Earth-like value γ (= 0.019) and comet-like γ (= 0.71) higher than the Earth's: The equilibrium values of T s and P CO 2 for M oc < 86M oc,⊕ are increased up to 326 K and 2.4 × 10 -1 bars for γ = 0.19 and S ⊙ , respectively (see Figs. 5a and5c).

In contrast, when the HP ice is present (M oc ≥ 86M oc,⊕ ), the negative feedback never works and, consequently, the snowball state is attained. This is because the seafloor temperature on the area under the sorbet region, where seafloor weathering works, is fixed at the melting temperature of ice and, thus, insensitive to the surface temperature. Although the reduction in effective weathering area reduces seafloor weathering rate (see Fig. 5b), it is found to have little impact on surface temperature because the seafloor weathering flux is significantly higher than the degassing flux.

Dependence on stellar insolation

We examine the dependence of planetary climate on stellar insolation. Since the runaway greenhouse limit, which controls the inner edge of the habitable zone, is only slightly higher than S ⊙ (e.g., 1.06S ⊙ Kopparapu et al. 2013), we show only the results for smaller stellar insolation of 0.9S ⊙ than the fiducial value of 1.0S ⊙ . As shown in Fig. 5c, stellar insolation affects CO 2 partial pressure both in the equilibrium and snowball states: the smaller the stellar insolation, the higher the CO 2 pressure is, as a whole: P CO 2 for S = 0.9S ⊙ is higher by a factor of ∼ 3 and by two orders of magnitude than that for S = 1.0S ⊙ in the equilibrium and snowball states, respectively. The other quantities are almost unaffected by stellar insolation. This is because the increase in CO 2 pressure compensates for the decline in stellar insolation so as not to change the surface temperature which controls weathering behavior and COM-HP in our climate model. Thus, variation in stellar insolation has little impact on planetary climate, provided the planet is located in the habitable zone.

Dependence on mean mantle heat flow

Next, we examine what impact the mean mantle heat flow q has on the surface and seafloor conditions. Figure 6 shows (a) the surface temperature, (b) effective weathering area, and (c) seafloor weathering flux for five different choices of q. The variation in mean mantle heat flow turns out to yield no change on the overall behavior, but quantitative modifications to the ocean mass dependence. First, as seen in Fig. 6a, when no HP ice is present, the larger the mean mantle heat flow, the higher the surface temperature is for a given ocean mass. The variation in q leads to a large difference in the surface temperature (up to 40 K). Also, the surface condition lapses into the snowball state at larger ocean mass for a larger q. That is because as q increases, the seafloor spreading rate A S increases (see Eq. [21]) and, thus, the degassing flux increases (see Eq. [18]), leading to higher surface temperature and larger critical ocean mass for forming the HP ice (COM-HP, see also Fig. 1). In Fig. 6c, the seafloor weathering flux is also found to increase by approximately an order of magnitude in response to the rise in the degassing flux.

As shown in Fig. 6b, the effective weathering area f oc starts to decrease from unity at larger ocean mass for larger mean mantle heat flow and is always unity until M oc = 200M oc,⊕ for q ≥ 100 mW m -2 . While the seafloor weathering flux changes with q (i.e., f oc ) in the case with HP ice, the reduction in f oc turns out to have only a small effect on the surface temperature, because of significantly high seafloor weathering rate for any value of q (see Fig. 6c).

Dependence on CO 2 /H 2 O ratio and seafloor weathering activation energy

As described in section 2.3, the carbon cycle depends on the total carbon inventory and seafloor weathering rate. The former may differ greatly from planet to planet, as suggested, for example, by a difference in the CO 2 /H 2 O molar ratio γ (Eq. [13]) between comets and the Earth. Also, the seafloor weathering rate is in general uncertain, mainly because the activation energy E a (Eq. [22]) is poorly determined observationally. Here we investigate the sensitivities of the planetary climate to γ and E a with focus on the critical ocean mass, beyond which the planetary climate is in the snowball state (hereafter, abbreviated to COM-SB and denoted by M cr(sb) oc

). In Fig. 7, we plot the relationships between M cr(sb) oc and total degassing flux, F D + F M , for various values of γ between 7.4 × 10 -3 and 2.1; both M cr(sb) oc and F D + F M are obtained from the carbon cycle calculations. Here we assume q = 80 mW m -2 . Fig. 7a shows the fiducial case with E a = 41 kJ mol -1 ; Fig. 7b shows cases with three different values of E a . For reference, in Fig. 7a, we show the relationships between F D + F M and not M cr(sb) oc but M oc for three different values of γ by dashed lines (the result for γ = 0.19 is also shown in Fig. 5b). In Fig. 7a, we can see that the total degassing flux increase almost linearly with γ for a given ocean mass. For γ < 7.4 × 10 -3 , the COM-SB is absent be- cause the snowball state is achieved in all the ocean mass range due to low degassing flux.

As shown in Fig. 7a, the total degassing flux has a peak at M cr(sb) oc = 128M oc,⊕ . For M cr(sb) oc ≤ 128M oc,⊕ , the COM-SB increases from 24 to 128 M oc,⊕ and the total degassing flux, which is determined by F SW (T floor , f oc ) = F SW (T mel floor , 1), increases from 8.1 × 10 11 to 6.8 × 10 14 mol yr -1 with increase in γ from 7.4×10 -3 to 2.1. Despite order-of-magnitude variation in F D + F M , the COM-SB varies moderately by a factor of ∼5 (see section 5.1 for an analytical interpretation). On the other hand, for M cr(sb) oc > 128M oc,⊕ , the total degassing flux is determined by the minimum weathering flux with the HP ice, namely F SW (T floor , f oc ) = F SW (T mel floor , f oc,min ) (see Fig. 4 for f oc,min ). Thus, the COM-SB increases from 128 to 200M oc,⊕ and the total degassing flux decreases from 6.8 × 10 14 to 2.3 × 10 14 mol yr -1 with decrease in γ from 2.1 to 5.2 × 10 -1 . In this diagram, equilibrium climates (F SW = F D + F M ) are achieved on the side above the solid line, whereas the planetary surface condition lapses into snowball states (F SW > F D + F M ), because of the presence of the sorbet region, on the side below the solid line. Note the curve of M cr(sb) oc is a bit inflected at M oc = 85M oc,⊕ because of a phase change of the HP ice.

In Fig. 7b, we show the impact of E a on M cr(sb) oc

. The curve for larger E a is found to be steeper. The three curves cross each other at F D + F M = 2.0 × 10 12 mol yr -1 , where the seafloor temperature T floor is equal to T 0 (=289 K), so that F SW is independent of E a (see Eq. [22]). Above the crossover point, the higher the activation energy, the smaller the COM-SB is; its dependence is opposite below the crossover point. Although being large relative to that on γ, the dependence of M cr(sb) oc on E a is at most linear. Thus, it would be fair to say that M cr(sb) oc is rather insensitive to E a . We further discuss the nature of the COM-SB analytically in section 5.1.

DISCUSSION

Critical ocean mass for snowball state

One of the most important findings in this study is that there is a critical ocean mass, beyond which an ocean terrestrial planet has an extremely cold climate (i.e., the snowball state). Furthermore, we have found that the COM-SB,

M cr(sb) oc

, falls into a relatively narrow range between 20 and 100 M oc,⊕ . Here we give an interpretation to the low sensitivity of M cr(sb) oc to the planetary mass M p , the total degassing flux F D + F M , and the activation energy of seafloor weathering E a , by deriving an approximate solution for M cr(sb) oc . This could help us obtain an integrated view of planetary climate on ocean planets under our idealized seafloor environments.

As demonstrated in section 4, the planetary climate lapses into the extremely cold one, when HP ice is formed on the seafloor. Then, the seafloor temperature T floor is fixed at the melting temperature T mel and thus determined uniquely by the seafloor pressure P floor . Since the ocean mass and depth are negligibly small relative to the planetary mass and radius, respectively, under hydrostatic equilibrium, the seafloor pressure is given approximately by

P floor ≈ GM p 4πR 4 p M oc = G ρp 3R p M oc , (23)
where R p and ρp are the planetary radius and mean density, respectively. For M oc = M cr(sb) oc

, P floor corresponds to the crossover pressure between the adiabat and the melting curve, both of which are independent of planetary mass. Thus, from Eq. ( 23 [22]). Panel (a) shows the fiducial case with E a = 41 kJ mol -1 . For reference, dashed lines represent the relationship between F D + F M and not the COM-SB, but just the ocean mass for three different values of γ when the planetary climate is in an equilibrium state. Panel (b) compares the results for three different values of E a . Dashed lines represent the analytical solutions of COM-SB (see Eq.

[25]). On the right side of the solid lines, the sorbet is present in the deep ocean. In these calculations, we have assumed the mean mantle heat flux q = 80 mW m -2 .

planetary mass; indeed, between M p = 1M ⊕ and 10M ⊕ , for example, M cr(sb) oc differs only by ∼12%. To derive the dependence of M cr(sb) oc on F D + F M and E a , we consider seafloor weathering. In the equilibrium state, since F SW (T floor ) = F D + F M , T floor is given as a function of Eq. [22]). Also, T floor = T mel , when M oc = M cr(sb) oc : From Eq. ( 5),

F D +F M (see
T mel ≃ c 1 + c 2 P floor , (24)
where c 1 = 236 K, c 2 = 6.09 × 10 -8 K Pa -1 . From Eqs. ( 22)-( 24), M cr(sb) oc is expressed as

M cr oc = 3R p c 2 G ρp ⎧ ⎨ ⎩ T 0 1 -RT 0 E a ln F D +F M f oc F * SW -c 1 ⎫ ⎬ ⎭ . (25)
This equation confirms that the COM-SB depends on the degassing flux only weakly. Also, since the denominator of the first term must be positive, the sensitivity of M cr(sb) oc to E a turns out to be small. In Fig. 7b, we plot the relationships between M cr(sb) oc and F D + F M calculated from Eq. ( 25), which is found to reproduce the numerical results well, except for the effect of phase change of the HP ice.

Effect of supply limit of cations

As shown in section 4, without any limit to seafloor weathering, the presence of the HP ice (exactly to say, the sorbet) always enhances seafloor weathering, resulting in extremely cold climates (i.e., the snowball states). In reality, however, the seafloor weathering rate is limited by the number of cations available in the oceanic crust. This is because seafloor weathering occurs through hydrothermal circulation The horizontal black solid line represents the supply limit above which seafloor weathering F SW is limited by the insufficient supply of cations (see Eq. [26]). The green solid line corresponds to the critical ocean mass (COM-SB), namely, the boundary between the equilibrium states and the extremely cold (snowball) states that would be achieved if no supply limit were assumed, same as that in Fig. 7 a. In this calculation, we have assumed the mean mantle heat flux q = 80 mW m -2 and activation energy E a = 41 kJ mol -1 .

in the oceanic crust and thus the amount of cations available depends on the depth of hydrothermal circulation. This limit to seafloor weathering rate, which we call the supply limit, F limit SW , can be given by (Sleep et al. 2001;Foley 2015)

F limit SW = x m c ρ rock d hy A S , (26)
where x is the number fraction of cations (Ca 2+ , Mg 2+ , and Fe 2+ ) in the oceanic crust, m c is the averaged molar mass of cation, ρ rock is the density of the oceanic crust, d hy is the depth at which hydrothermal carbonation occurs, and A S is the seafloor spreading rate. In the present Earth condition, where x = 0.3, m c = 55 g mol -1 , ρ rock = 2800 kg m -3 and d hy = 500 m (Sleep et al. 2001), F limit SW = 7.6 × 10 6 A S . If the total degassing flux is higher than the supply limit, the atmospheric CO 2 continues to increase with age. Qualitatively, the more the atmospheric CO 2 , the higher the surface temperature is. While climate sensitivity to the amount of CO 2 is unclear for high P CO 2 because of poor understanding of radiative forcing of water vapor for hot atmospheres, recent 1-D radiative-convective calculations show that T s > 350 K for P CO 2 > several bars, provided the stellar insolation is equal to that for the present Earth (Wordsworth & Pierrehumbert 2013;Ramirez et al. 2014).

Figure 8 is the climate diagram for q = 80 mW m -2 and E a = 41 kJ mol -1 , where we indicate three different climate regimes, which include the equilibrium climates, the extremely cold climates (or the snowball states), and the extremely hot climates. The extremely hot climate is a state such that F SW < F D + F M , because of supply limit so that CO 2 accumulates in the atmosphere. The supply limit (horizontal black solid line) is calculated from Eq. ( 26). The boundary between the equilibrium-climate and cold-climate regimes (green line) corresponds to the COM-SB shown in Fig. 7a. Of importance here is that the total degassing flux at the COM-SB is always higher than the supply limit for M oc > 70M oc,⊕ . Thus, for M oc > 70M oc,⊕ , the planet has no equilibrium climate (i.e., extremely hot or cold climate) because of the enhanced seafloor weathering and the supply limit.

Here we give a brief discussion about the uncertainty in the supply limit. Although the mean mantle heat flow q, which determines the seafloor spreading rate A S and thus the supply limit F limit SW , decreases with age during planetary evolution, its decrement on a timescale of billion years is known to be similar to the mean mantle flow for Earth-like planets with age of several billion years (McGovern & Schubert 1989). Also, we have adopted the solidus temperature of rock for T sol in Eq. ( 11), instead of the potential temperature of the mantle, which leads to overestimating the supply limit approximately by a factor of 2 in the case of T sol = 2000 K, which corresponds to the potential temperature at hot initial states (e.g., [START_REF] Tajika | [END_REF]. In addition, in the equilibrium states, the effects of variation in seafloor spreading rate are canceled out, because both of the supply limit and degassing flux have a linear dependence on seafloor spreading rate (Eqs. [18] and [26]). Thus, the uncertainty in q has a small influence on the climate diagram for ocean planets.

The hydrothermal carbonation depth d hy would depend on ocean mass. Some experiments suggest that the hydrothermal carbonation depth decreases with increasing seafloor pressure because thermal cracking becomes weaker (Vance et al. 2007). Thus, the supply limit is expected to decrease with ocean mass, which would extend the domain of the extremely hot climate in Fig. 8.

In conclusion, the enhanced seafloor weathering due to the formation of the sorbet region and the supply limit narrow the range of ocean mass of terrestrial planets with the equilibrium climates. This implies the difficulty of clement climates, like the present Earth, on ocean planets with plenty of water.

Caveats

Ocean Layer Model

Here we discuss the validity of our assumptions regarding the ocean layer, which include: (1) No boundary layer exists at the top of the HP ice layer; (2) The carbon partitioning between the atmosphere and ocean is always in equilibrium and the CO 2 content is constant through the sorbet region; and (3) The heat transport occurs in the vertically one dimension.

(1) Regarding convective transport in the HP ice, we have considered the presence of a thermal boundary layer at the bottom, but not at the top. To evaluate the effect of the top boundary layer (TBL) on the effective weathering area f oc , we have calculated f oc in the same settings as in Fu et al. (2010), who considered TBL in addition to a bottom boundary layer (BBL). Then, we have found that TBL leads to reducing the effective weathering area in the low surface temperature domain for a given ocean mass (e.g., T s 390 K for M oc = 200M oc,⊕ , see also Fig. 3 f ). This is because BBL is cooler and thicker without TBL than in the case with TBL. As discussed in section 3.2, in this domain, the reduced thickness of BBL, δ , increases the critical heat flow (Eq. [10]) and, thus, reduces the effective weathering area (Eq. [12]). However, we have also found that the presence of TBL brings about little change in the maximum of q cr at a given ocean mass (Fig. 4). This is because the same temperature gradient in BBL is achieved by a change in T s , given that TBL is assumed to follow the melting line of ice (see Fig. 2 in Fu et al. 2010). Hence, the climate diagram for ocean planets is almost unaffected by the presence of TBL.

(2) We have assumed that the CO 2 circulation in the sorbet region occurs efficiently enough that carbon partitioning between the atmosphere and ocean remains in equilibrium. However, the CO 2 circulation (not the seafloor weathering) limits the consumption of atmospheric CO 2 , if being slower than the response of the carbon budget in the oceanatmosphere system. The latter is controlled by regassing in the environments of interest in this study, although depending on degassing, in general [START_REF] Tajika | [END_REF]; thus, its timescale is ∼ τ/β .

The CO 2 circulation occurs in the following way: Aqueous CO 2 converts to CO 2 ice quickly within the HP ice (Bollengier et al. 2013) and, then, the CO 2 ice moves together with the HP ice. In the HP ice layer, since the upward sorbet flow transports mass (and heat), the HP ice sinks accordingly for mass conservation. Below we estimate the sinking speed of the HP ice and the overturn timescale of the HP ice layer from energy balance and mass conservation.

When heat is transported by thermal diffusion and sorbet flow, the energy balance is expressed as

Q = -k dT dr + χρ l (L +C P ∆T HP )w l , ( 27
)
where Q is the heat flux from the oceanic crust, χ is the melt fraction, L is the latent heat of HP ice, and ρ l , C P , and w l are the density, specific heat, and flow speed of liquid water in the sorbet, respectively. The first and second terms on the right side represent the thermal conduction and melt advection, respectively. Note that we have assumed that permeable flow of liquid dominates the sorbet flow and, namely, neglected upwelling solid flow, which results in underestimating the sinking flux of the HP ice. From mass conservation and Eq. ( 27), the sinking speed of the HP ice, w HP , against the upwelling sorbet flow is given by

w HP = χρ l (1 -χ)ρ HP w l = Q + k dT dr (L +C P ∆T HP )(1 -χ)ρ HP , (28)
where ρ HP is the density of the HP ice. The thermal conduc- tion flux along the melting curve is -kdT /dr ≈ 10 mW m 2 (see Fig. 3b). The heat flux from the oceanic crust of 80 mW m -2 being added, Q + kdT /dr ≈ 70 mW m -2 . The material properties of liquid water and HP ice are ρ HP = 1400 kg m -3 , L = 4.2 × 10 5 J kg -1 (at 300 K from Dunaeva et al. 2010), and C P = 4.1 × 10 3 J kg -1 K -1 (at 300 K from Waite et al. 2007).

For terrestrial sea ice, permeability decreases abruptly for melt fraction below χ = 5 % (Golden et al. 1998). Although not known well for the HP ice, we assume that the HP ice behaves in a similar way and use χ = 5 %. According to our calculation results, ∆T HP = 78 K and D = 99 km for M oc = 100 M oc,⊕ and T s = 300 K. Then, the overturn timescale of the HP ice (≡ D/w HP ) comes out to be 44 Myr. Even for the M oc = 200M oc,⊕ , D/w HP = 203 Myr. On the other hand, τ/β ∼ 250 Myr for the present Earth's condition, using the value of τ for the present Earth (∼ 100 Myr) [START_REF] Turcotte | Geodynamics -2nd Edition[END_REF]. Thus, the CO 2 circulation occurs faster than the response of the carbon budget in the atmosphere-ocean system.

The above estimate may remain to be refined. For example, using the relation between the Nusselt number (Nu) and the Rayleigh number, Nu ∝ Ra 1/3 , [START_REF] Turcotte | Geodynamics -2nd Edition[END_REF], the semi-infinite half-space cooling model shows that the residence timescale τ is ∝ η 2/3 man , where η man is a viscosity of the mantle material, and, thus, depends strongly on mantle temperature and water content in the mantle. Indeed, the residence timescale is thought to have varied by an order of magnitude during the thermal evolution of the Earth [START_REF] Tajika | [END_REF]. Also, seafloor weathering would be limited, if the planet has a thick HP ice and vigorous convective mantle. However, it is emphasized here that given a weak dependence of the COM-SB even a sluggish circulation of CO 2 in the HP ice could yield no significant change in COM-SB.

(3) We have assumed a vertically one-dimensional structure of the ocean and thus considered only vertical heat transport. In reality, the thermal structure of the ocean is more complicated because of convective patterns and inhomogeneous phase change. First, since the distance between the mid-ocean ridge and trench ( 10000 km) is much larger than the thickness of the HP ice (100 km), with which the size of convective cells is comparable (e.g., [START_REF] Turcotte | Geodynamics -2nd Edition[END_REF], detailed convective patterns matter little for the overall heat transfer in the HP ice layer. In addition, hydrodynamical simulations show a heatpipe structure of the HP ice layer for high heat fluxes from the oceanic crust, which means phase change rarely occurs vertically throughout the ocean (Choblet et al. 2017;Kalousová et al. 2018).

Atmospheric model

Our climate modeling has demonstrated that the runaway cooling due to atmospheric CO 2 drawdown generally occurs on ocean planets with plate tectonics, provided M oc > M cr(sb) oc . Although we assume the runaway cooling ends up with the snowball state with T s < 273 K, it is to be examined more carefully whether the global snowball state is achieved or not. Here we discuss some uncertainties of the atmospheric model.

We have assumed a constant surface albedo of 0.32. Planetary albedo depends on cloud radiative forcing that generally depends on T s (e.g., Wolf & Toon 2013, 2015). As far as partially ice-covered planets are concerned, simulations based on 3-D general circulation models (GCMs) for the Archean Earth (Wolf & Toon 2013) and ocean planets (Charnay et al. 2017) demonstrate that planetary albedo increases rapidly with decreasing T s due to the ice-albedo feedback, although the contribution of clouds to planetary albedo declines. This means that the assumed surface albedo of 0.32 is an underestimate for the planetary albedo for cold climates of interest in this study and thereby leads to overestimating T s . This indicates that the runaway cooling results in the snowball state, even if the planet receives high stellar insolation comparable to the present Earth.

Also, we have assumed that the distribution of relative humidity in the atmosphere is the same as that in the present Earth's atmosphere. Wordsworth & Pierrehumbert (2013) and Ramirez et al. (2014) found multiple equilibrium solutions for a CO 2 -free, almost water-saturated atmosphere, including a hot solution with the surface temperature of ∼ 500 K, even if stellar insolation is comparable to the present Earth's stellar insolation. This suggests that the snowball state is not always achieved. However, GCM simulations show that atmospheric circulation leads to precipitation and thereby to removing water vapor from the atmosphere, namely, making unsaturated regions even if stellar insolation is close to the runaway greenhouse limit (Wolf & Toon 2015). This suggests that such a hot state would be unlikely to occur, although more work is needed to confirm so.

Carbon cycle model

In this study we have ignored the situation where the atmosphere is so cold that the surface of the ocean is frozen and, instead, have stopped calculations once the surface temperature reaches 273 K. Here we discuss the effect of surface ice on the carbon cycle and warming process in the snowball state. When the surface ice is convectively stable, which is appropriate for Earth-like high heat fluxes and moderately low surface temperatures (Fu et al. 2010), molecular diffusion in the surface ice would control the exchange of CO 2 between the atmosphere and ocean. Performing molecular dynamics simulations, Ikeda-Fukazawa et al. (2004) estimated that CO 2 molecular diffusion coefficient in H 2 O ice is ∼ 10 -10 m 2 s -1 at 270 K. For the thickness of surface ice of 1 km, for example, the diffusion timescale is on the order of Gyr. This means that even for a degassing flux higher than the critical value shown in Fig. 7, CO 2 accumulation in the atmosphere proceeds too slowly for the climate to escape from the snowball state. This indicates that the snowball state we have found is maintained on a timescale of Gyr. We have to keep in mind, however, that it still remains a matter of debate how past Earth escaped from the snowball state.

On the other hand, the seafloor weathering is thought to be insensitive to the existence of surface ice, as follows. The seafloor temperature is fixed at the melting temperature of the ice in the snowball solutions. Since surface ice has a steep conductive temperature gradient and thus the thickness is small relative to the whole ocean, the T -P structure below the surface ice is rather insensitive to the surface temperature, thereby having little effect on the seafloor pressure and temperature. Thus, the seafloor weathering flux is proportional to the effective weathering area. The latter increases with decreasing the surface temperature in the low surface temperature regime shown in Figure 3 f . Thus, beyond the critical ocean mass, the seafloor weathering would be higher than the degassing flux, even if the surface ice is formed. Thus, once being achieved by the runaway cooling, the snowball state is maintained.

Exoplanet

Finally, we discuss an application of our findings to terrestrial exoplanets. Although we have no enough knowledge of the degassing flux of exoplanets, which depends on several uncertain factors such as planetary carbon budget, thermal structure of planetary interior, and ocean mass, we have found that the COM-SB is less sensitive to the degassing flux. As shown in Fig. 7a, terrestrial exoplanets with oceans of more than several tens of M oc,⊕ in the habitable zone have extremely cold climates. Cold climates are also suggested for Earth-like planets with low degassing flux in the habitable zone (e.g., Kadoya & Tajika 2014). Thus, terrestrial planets with CO 2 -poor cold climates would not be uncommon in the habitable zone around Sun-like stars, provided plate tectonics is common for those planets.

Recently, habitability for planets around ultra cool stars (e.g., Proxima Centauri and TRAPPIST-1) are actively debated (e.g., Ribas et al. 2016;Turbet et al. 2018;Valencia et al. 2018). Since the snowline is located near the habitable zone and ice-rich planets readily migrate from beyond, ocean planets would be abundant in the habitable zone around cool stars (e.g., Tian & Ida 2015). Planets around cool stars are synchronously rotating, which results in a large difference in surface temperature between the day and night sides (e.g., Pierrehumbert 2011). This might result in different phase structure and flow pattern in the ocean layer from our situation and CO 2 condenses on the night side (Turbet et al. 2018). In this case, the HP ice would be easily formed on the cool night side.

However, provided all our assumptions are valid also for synchronously rotating planets and the dayside and nightside have the same thickness of the ocean layer, the weathering flux on the dayside is always higher than that on the nightside, because of high surface temperature due to the concentration of all the stellar insolation. Thus, an equilib-rium climate could be achieved on the dayside, although the nightside is extremely cold. Then, the COM-SB for such a planet can be defined in the same way as we have done above and its value is equivalent to the estimate given in the previous sections. This implies a low probability of exoplanets with temperate climates in the habitable zone also around cool stars. Note that even if they have massive oceans with a mass larger than the COM-SB, synchronously rotating planets never become snowballs, because the local climate around the substellar point could be always temperate (Checlair et al. 2017).

SUMMARY AND CONCLUSION

The Earth's climate is stabilized by temperature-dependent, efficient continental weathering. Beyond the solar system, however, there must be continent-free terrestrial planets covered with global oceans (called ocean planets). Only with inefficient seafloor weathering, the Earth's climate would be much warmer. Furthermore, previous studies suggest that ocean planets have extremely hot climates, if they have massive oceans of 20 to ∼ 100M oc,⊕ , because the HP ice present in the deep ocean completely prevents chemical weathering on the oceanic crust (Alibert 2014;Kitzmann et al. 2015). However, those studies oversimplify the heat transfer in the HP-ice layer and ignore horizontal variation from heat flow from the oceanic crust. Thus, in this study, we have revisited the climate of ocean planets with plate tectonics in the habitable zone, by incorporating the effects of the liquid/solid coexistence region (called the sorbet region) near the midocean ridge in the carbon cycle (Fig. 2). The main findings of this study are summarized as follows.

Our seafloor environment model without the effect of the carbon cycle (i.e., fixed surface temperature) has shown that even if pressures in the deep ocean are high enough for HP ice to form, heat flux from the crust is too high to be transferred by solid convection, making the HP ice molten and forming a sorbet region, at least, near the midocean ridge (section 3). Although reduced with increasing ocean mass or decreasing mean mantle heat flow, the effective weathering area never becomes zero for M oc ≤ 200M oc,⊕ . This means that seafloor weathering remains possible and subsequent material circulation (e.g., carbon cycle) will sufficiently occur through the sorbet region.

Modeling the carbon cycle with the effect of seafloor weathering under the sorbet region, we have found that the climate on the ocean planet is destabilized and lapses into a CO 2 poor, extremely cold state, which is called the snowball state (section 4). Such destabilization is triggered because seafloor temperature is fixed at the melting temperature of the HP ice and, thus, a high seafloor weathering flux is kept regardless of surface temperature, unlike continental weathering which is dependent on surface temperature. This indicates the existence of a critical ocean mass, beyond which an ocean planet no longer maintains a temperate climate. We have demonstrated that the critical ocean mass is less sensitive to planetary mass, degassing flux, and the detailed dependence of seafloor weather flux on seafloor temperature (i.e., the activation energy E a ), and is several tens of M oc,⊕ . Also, because of the supply limit of cations, seafloor weathering is ineffective in compensating massive degassing, not achieving equilibrium climates, but yielding extremely hot ones.

As demonstrated in this paper, thermal and chemical interaction between the ocean and rocky interior significantly alters the planetary climate of ocean planets even in the habitable zone. We have found that temperate equilibrium climates are achieved in limited ranges of ocean mass and degassing flux. This suggests that a certain proportion of terrestrial exoplanets in the habitable zone could be frozen ocean planets, provided they are Earth-like ones with plate tectonics. In any case, our findings indicate that ocean mass has a crucial role in the planetary climate of terrestrial planets with a massive ocean. While the characterization of terrestrial exoplanets will be performed for detecting habitable planets in the next decade, we should discuss their climates carefully because those exoplanets would be diverse in surface water amount.

We have assumed the Henry's law is valid. We use the values of K 0 from Weiss (1974), K 1 and K 2 from Millero et al. (2006), and K w from Millero (1995), which were obtained experimentally for temperature of 300 K and salinity of 35 % . We use those values at all temperatures, because those equilibrium constants have not been experimentally measured for the temperature range of interest in the study.

Since Equations (B5)-(B8) contain six unknowns, we need at least two additional equations. One is the equation of charge conservation:

[OH -] + [HCO - 3 ] + 2[CO 2- 3 ] = [H + ] + [M + ], (B9)
where [M + ] represents the concentration of all the cations in the ocean. In this study, we use the average value of [M + ] measured in the present Earth ocean (i.e., [M + ] = 2.2 × 10 -3 mol L -1 ). Even a constant value of [M + ] has little influence on conclusions (see also section 2.3). Also, the total number of carbon C atm +C oc , which is determined from the carbon cycle model (section 2.3), must be conserved in the atmosphere-ocean system. Since the solubility of CO 2 in liquid water increases rapidly with pressure (Duan & Sun 2003) and mixing occurs in the ocean on a timescale much shorter than that of interest in this study, we assume C oc to be equal to the surface concentration of CO 2 . Thus, the total number of carbon is expressed as

C atm +C oc ≈ 4πR 2 p P CO 2 m CO 2 g s + ([H 2 CO 3 ] + [HCO - 3 ] + [CO 2- 3 ])V oc , (
B10) where the first term on the right-hand side corresponds to C atm , and m CO 2 is the molecular weight of CO 2 (= 44 g mol -1 ) and g s is the surface gravity. We obtain R p , g s , and V oc for given T s and M oc from the internal structure model in section 2.1. For C atm , we assume that the molecular weight of the atmospheric gas is equal to the molecular weight of CO 2 , which overestimates P CO 2 . However, the approximation has no influence on the overall results of the study.

Finally, solving Eqs. (B5)-(B10), we determine P CO 2 and the mole fractions of ions.

Figure 2 .

 2 Figure2. Seafloor environment model-(a) schematic illustration of the seafloor environment and (b) qualitative temperature profiles in the infinitesimally thin layer on the side of the high-pressure ice at the boundary between the "sorbet" and "high-pressure ice" regions. In panel (a), arrows represent the direction and dominant mechanisms of heat transport. In the sorbet region, ice and liquid coexist and thus the temperature is fixed at the melting point of H 2 O. In panel (b), Ra and Ra cr represent the Rayleigh number and the critical Rayleigh number, respectively. The red and orange solid lines represent thermal structure when the high-pressure ice layer is convective and conductive, respectively. The blue and black solid lines represent the melting line of H 2 O and adiabatic thermal structure of the liquid ocean, respectively.

Figure 3 .Figure 4 .

 34 Figure3. Formation of high-pressure (HP) ice and its impacts on the seafloor condition. Thickness of the HP ice layer (panels a and d), critical heat flow (panels b and e), and effective weathering area (panels c and f ) are shown as a function of ocean mass M oc for surface temperature T s = 300 K (top) and as a function of T s for M oc = 200M oc,⊕ (bottom). M oc,⊕ represents the present Earth's ocean mass. Note that we have not used the carbon cycle model for determining T s here, but performed calculations for given values of T s , instead.

Figure 5 .

 5 Figure 5. Surface and seafloor conditions obtained from the carbon cycle model: (a) Surface temperature (red sold line) and seafloor temperature (red dashed line), (b) seafloor weathering flux (blue solid line) and effective weathering area (blue dashed line), and (c) partial pressure of CO 2 are shown as a function of ocean mass in the unit of the Earth's ocean mass M oc,⊕ . Shaded is the range for the snowball state. In this calculation, the mean mantle heat flow is assumed to be 80 mW m -2 . The symbol S ⊙ represents the solar insolation received by the present Earth.

Figure 6 .

 6 Figure 6. Surface and seafloor conditions obtained from the carbon cycle model for five different values of the mean mantle heat flow. (a) Surface temperature, (b) effective weathering area and (c) seafloor weathering flux are shown as functions of ocean mass.

,Figure 7 .

 7 Figure7. Relationship between the critical ocean mass for the snowball state (COM-SB) and the total-degassing flux F D + F M (solid lines) for different choices of the CO 2 /H 2 O molar ratio in the source of volatiles, γ (see Eq. [13]), and seafloor weathering activation energy, E a (seeEq. [22]). Panel (a) shows the fiducial case with E a = 41 kJ mol -1 . For reference, dashed lines represent the relationship between F D + F M and not the COM-SB, but just the ocean mass for three different values of γ when the planetary climate is in an equilibrium state. Panel (b) compares the results for three different values of E a . Dashed lines represent the analytical solutions of COM-SB (seeEq. [25]). On the right side of the solid lines, the sorbet is present in the deep ocean. In these calculations, we have assumed the mean mantle heat flux q = 80 mW m -2 .

FFigure 8 .

 8 Figure 8. Climate diagram that shows different climate regimes on the plane of the total degassing flux (F D + F M ) vs. ocean mass.The horizontal black solid line represents the supply limit above which seafloor weathering F SW is limited by the insufficient supply of cations (seeEq. [26]). The green solid line corresponds to the critical ocean mass (COM-SB), namely, the boundary between the equilibrium states and the extremely cold (snowball) states that would be achieved if no supply limit were assumed, same as that in Fig.7a. In this calculation, we have assumed the mean mantle heat flux q = 80 mW m -2 and activation energy E a = 41 kJ mol -1 .

Table 1 .

 1 Coefficients

		for ice melting curve given by Eq. (5) from
	Dunaeva et al. (2010)
	Ice phase	a 1	a 2	a 3	a 4	a 5
	Ice VI	4.2804	-0.0013	21.8756	1.0018	1.0785
	Ice VII	-1355.42	0.0018	167.0609	-0.6633	0

Table 2 .

 2 Variables and their values.

	Parameter	Symbol	Value
	Ocean mass	M oc	1-200 M oc,⊕
	Mean mantle heat flow	q	40, 60, 80, 100, 120 mW m -2
	Activation energy of seafloor weathering	E a	30, 41, 92 kJ mol -1
	CO 2 /H 2 O molar ratio	γ	1.0 × 10 -3 -10

Table 3 .

 3 Parameters and their values.

	Parameter	Symbol	Value
	Earth ocean mass Molar quantity of H 2 O in the Earth ocean mass Thermal conductivity of the HP ice	M oc,⊕ n oc,⊕ k	1.37 × 10 21 kg 7.6 × 10 22 mol 3.8 W m -1 K -1
	Constant for the viscosity of the HP ice Characteristic shear stress of the HP ice Activation energy for the viscosity of the HP ice	B ζ E *	6.7 × 10 19 Pa 4.5 s 2.0 × 10 6 Pa 110 kJ mol -1
	Activation volume for the viscosity of the HP ice Critical Rayleigh number	V * Ra cr	1.1 × 10 -5 m 3 mol -1 2000
	Thermal conductivity of the oceanic crust	k rock	3.3 W m 2 s -1
	Thermal diffusivity of the oceanic crust Present Earth's seafloor weathering rate Reference seafloor temperature	κ rock F * SW T 0	1.0 × 10 -6 m 2 s -1 2.0 × 10 12 mol yr -1 289 K
	Area of the oceanic floor Regassing ratio	A 0 β	3.1 × 10 14 m 2 0.4

MNRAS 000, 1-19 (2019)

https://github.com/VirtualPlanetaryLaboratory/atmos MNRAS 000, 1-19 (2019)

ACKNOWLEDGEMENTS

The authors thank the anonymous reviewer for thoughtful comments that greatly improved the manuscript. This work was supported by JSPS KAKENHI No. JP18H05439, JP23103003, and JP17H06104, ABC-NINS No. AB301002, JSPS Core-to-core Program "International Network of Planetary Sciences", and the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 679030/WHIPLASH).

APPENDIX A: INTERNAL STRUCTURE MODEL

We develop a radially one-dimensional, hydrostatic internal structure model in section 2.1, based on Valencia et al. (2007a). We consider a differentiated solid rock-metal body of 1 Earth mass covered with various amounts of H 2 O. Note that the planetary mass is the sum of the rock-metal body and ocean masses (i.e., 1M ⊕ + M oc ). We assume (1) the mass ratio of iron core to silicate-mantle is 7 : 3, (2) the mass ratio of the inner to outer core is the same as that of the Earth (35 : 65), (3) phase transitions occur at the same pressures as in the Earth's interior, and (4) thermal expansion of the mantle and core never occurs. These assumptions have little influence on the surface gravity, which affects the structure of the ocean layer and thus on our conclusion in this study.

The equations of state (EOSs) and parameter values that we adopt are summarized in Table A1. The temperature effect on density of the HP ice follows expressions from Bezacier et al. (2014) for ice VI and Fei et al. (1993) for ice VII. The thermal capacity of liquid water is taken from Waite et al. (2007). The phase transition from water to HP ice occurs where the adiabat crosses the phase boundaries in the P-T plane given by Eq. 5.

We assume the thermal structure of the iron core and the rocky mantle are the same as the Earth. The phase transitions and the pressures at the transitions are summarized in Table A2.

APPENDIX B: PARTITIONING OF CO 2 BETWEEN ATMOSPHERE AND OCEAN

The partial pressure of CO 2 , P CO 2 , depends on the carbon budget of the surface reservoirs (C atm + C oc ), oceanic pH and ocean volume, V oc . Here we describe the calculation method of carbon partitioning between the atmosphere and ocean, which is almost the same as that described in [START_REF] Tajika | The evolution of the terrestrial environment[END_REF] and Kitzmann et al. (2015).

Chemical equilibrium among CO 2 (g), carbonic acid (H 2 CO 3 ), carbonate ion (CO 2- 3 ), and bicarbonate ion (HCO - 3 ) is determined by the following reactions

with equilibrium constants

)

)

)

APPENDIX C: DEPENDENCE OF DEGASSING COEFFICIENT ON SEAFLOOR PRESSURE

Here we introduce the dependence of degassing coefficient, K D , on seafloor pressure. According to [START_REF] Tajika | [END_REF], the degassing coefficient is given by

where f CO 2 is the degassing fraction, which is defined as the molar fraction of the CO 2 degased from the upwelling magma at the ridge, d m is the degassing depth, which is defined as the melt generation depth of mantle, and V man is the volume of the mantle. For V man and d m , we use the values for the present Earth, namely V man = 8.0 × 10 20 m 3 and d m = 40 km [START_REF] Tajika | [END_REF].

The degassing fraction f CO 2 depends on the ocean mass, because of pressure dependence of CO 2 solubility into magma (Kite et al. 2009). In this study, we take it into account, following [START_REF] Tajika | [END_REF], who considered the solubility equilibrium of CO 2 with solid/liquid silicate. We incorporate the pressure dependence on the solubility of CO 2 into silicate melts, K G/L , and the molar volume of CO 2 , V M , in calculating f CO 2 as [START_REF] Tajika | [END_REF]

where f melt is the melt fraction, K L/S is the partition coefficient of CO 2 between solid and liquid, and w G /w L represents the mass ratio of CO 2 partitioned into the gas phase to that into the liquid phase (liquid phase meaning CO 2 dissolved in melt); w G /w L is defined as

where φ is the vesicularity of melt, ρ rock is the density of oceanic crust, and n CO 2 is the molar concentration of CO 2 gas in the vesicles. We adopt values of f melt , K L/S , φ , and n CO 2 from Tajika & Matsui (1992) Recent molecular dynamics simulations (Guillot & Sator 2011) predict higher solubility of CO 2 for > 2 GPa than that obtained according to Henry's law. Those simulations found an almost linear dependence on pressure and weakly correlation with temperature. We have estimated the relationship between K G/L and P based on tabular data for T = 1673 K and MORB composition presented in Guillot & Sator (2011):

Here P is the pressure in the unit of GPa. We evaluate V M at the seafloor pressure using the EOS based on molecular dynamics simulations (Duan & Zhang 2006), which is of wide application (i.e., up to 10 GPa and 2573.15 K). The temperature in the EOS corresponds to the solidus of anhydrous peridotite at the seafloor pressure, which is parameterized by Hirschmann et al. (2009).

Figure C1 shows the degassing fraction f CO 2 as a function of ocean mass M oc for T s = 300 K. f CO 2 is found to decrease with M oc , because the solubility of CO 2 increases with pressure. f CO 2 varies from 0.23 to 0.1 between 1 to 200 M oc,⊕ . At M oc = 77 M oc,⊕ , the slope of f CO 2 changes because seafloor pressure becomes higher than 2 GPa. Kite et al. (2009) proposed that degassing could be completely suppressed (i.e., f CO 2 = 0) for a 100 km ocean (roughly 40 M oc,⊕ in our model) because of higher solubility of CO 2 . In contrast, Fig. C1 indicates that degassing also occurs for larger M oc . Higher solubility would lead to no partitioning into the gas phase (w G /w L → 0). In this case, degassing fraction would be determined by two-phase partitioning between solid and liquid and consequently f CO 2 becomes 0.096. Therefore, our model results in degassing that mainly occurs as the liquid phase at high pressures. Also, in Fig. C1, f CO 2 is estimated to be 0.23 for M oc = 1 M oc,⊕ corresponding to seafloor pressure of 27 MPa, which is relatively smaller than the value (0.32) estimated according to the Henry's law, by [START_REF] Tajika | [END_REF]. This difference is due to higher solubility (216 ppm at 27 MPa) than that (100 ppm) of [START_REF] Tajika | [END_REF]. Note that low-pressure experiments suggest higher solubility than our model (Jendrzejewski et al. 1997). In any case, because a pressure range much higher than 27 MPa is of special interest in this study, we neglect this difference.

This paper has been typeset from a T E X/L A T E X file prepared by the author.