Thermocline Depth on Water-rich Exoplanets

Water-rich exoplanet (alternatively called as ocean planet) is a type of proposed terrestrial planet that is covered by a global ocean and without any exposed continent at the surface (Kuchner, 2003; Léger et al., 2004; Selsis et al., 2007; Sotin et al., 2007; Marcus et al., 2010; Nettelmann et al., 2011; Zeng & Sasselov, 2014; Thomas & Madhusudhan, 2016; Auclair-Desrotour et al., 2018; Kite & Ford, 2018). The ocean can reach tens of to hundreds of kilometers, which is much deeper than the global-mean depth of the ocean on Earth, $\sim$3 km (Stewart, 2008). Due to the deep ocean, the pressure at the bottom of the ocean can be extremely high and it results in the occurrence of high-pressure ice at the sea floor. The existence of water-rich exoplanets and their possible bulk compositions have not yet been demonstrated directly by observations, but implied by the mean density derived from observed masses and radii. Their rather larger radii and lower densities compared to other rocky planets (such as Earth) suggest that they may be volatile-rich or water-rich outside the rocky core (Zeng et al., 2019). Water-rich exoplanets may be common in the galaxy, especially around M-type stars (Mulders et al., 2015; Brugger et al., 2016; Alibert & Benz, 2017).

Carbonate-silicate cycle plays an important role in regulating planetary climate and can strongly affect planetary habitability (Walker et al., 1981; Pierrehumbert, 2010). On water-rich exoplanets, due to the lack of exposed continents and continental weathering, the carbonate-silicate cycle has a weaker dependence on atmospheric CO₂, which might make the climate being better at resisting changes in external forcings (Hayworth & Foley, 2020). Krissansen-Totton et al. (2021) suggested that the high pressure at the sea floor is unfavorable for magmatic outgassing and seafloor weathering. Considering the temperature dependence of seafloor weathering might improve the habitability of water-rich exoplanets (Abbot et al., 2012). Moreover, the carbon partition between the atmosphere and ocean might play a negative role in stabilizing the climate on water-rich exoplanets, due to the decreased capability of the ocean to absorb carbon dioxide with temperature increasing (Kitzmann et al., 2015).

Ocean circulation also can affect the habitability of water-rich exoplanets. Different from other terrestrial exoplanets with shallow oceans like Earth, there is no subtropical or subpolar gyres, boundary currents, and coastal upwelling motions on water-rich exoplanets due to the lack of exposed continents. In particular, coastal upwelling plays a significant role in the upward transport of nutrients from the deep ocean to the upper layers on Earth, which can directly affect the biological activity in the ocean (Hutchings et al., 1995). Olson et al. (2020) demonstrated the importance of considering the wind-driven upwelling for the nutrient transport and biological activity on exoplanets with exposed continents and land-sea contrasts. However, there is no coastal upwelling on water-rich exoplanets. The meridional overturning circulation can also transport heat poleward, regulate the mean climate and thereby play a significant role in affecting the planetary habitability. Vallis & Farneti (2009) suggested that the effect of different ocean depths on the overturning circulation and the meridional ocean heat transport is limited, as long as the ocean depth is deeper than the main thermocline (a layer where the vertical gradient of temperature is large and is not subject to seasonal variability; Stewart 2008). However, the maximum ocean depth they tested is 4 km, simulations with deeper oceans have not been investigated.

The thermocline can strongly influence buoyancy, circulation, and the exchange of oxygen, carbon, heat and other nutrients between the upper layer and the lower layer (Zelle et al., 2004; Cantin et al., 2011; Giling et al., 2017). The thermocline is a layer of water where the temperature decreases with depth more rapidly and the vertical temperature gradient is larger than it does in the layers above or below (Pedlosky, 2006; Stewart, 2008; Fiedler, 2010). The thermocline on Earth can be divided into two kinds, one is main thermocline, the other is seasonal thermocline that varies with the seasons (Stewart, 2008). The thermocline separates the warm mixed layer (turbulent and properties uniformly distributed vertically) from the cold deep water, and the mixing between the upper and lower layers is inhibited by the large gradient. Due to the lack of inter-region mixing, oxygen content below the thermocline rapidly depletes with depth increasing, since organisms utilize it and there is no source below the thermocline where there is no sunlight (Thistle, 2003). Given that density is partly determined by temperature, and the vertical density gradient within the thermocline is also larger than the upper mixed layer and the lower layer.

Here we firstly explore the depth of the thermocline under different ocean depths. Then, the sensitivity of the thermocline depth in wind-dominated system and mixing-dominated system to various planetary and oceanic parameters is investigated. Section 2 presents our model and experimental designs. Section 3 presents our results, including the wind-driven circulation (Section 3.1) and the thermal-driven circulation (Section 3.2). We summarize and discuss the results in Section 4.

The global ocean model used in this study is the Massachusetts Institute of Technology general circulation model (MITgcm; Marshall et al. 1997a, b). By default, we simulate one planet with a global ocean that has a uniform depth of 40 km and has no any continental barrier in the zonal direction. In meridional direction, the ocean covers from 80^∘S to 80^∘N; higher latitudes are not simulated and solid walls are placed at the southern and northern boundaries of our domain in order to avoid grid cell convergence at the polar points. The primitive equations in spherical coordinates we used are written as:

where $\vec{u}$ is the horizontal flow vector, and $u$ and $v$ are its zonal and meridional components, respectively; $w$ is the vertical current speed; $\theta$ and $S$ are the potential temperature and salinity, respectively; the pressure field is splitted into a reference function of height $p_{0}(z)$ and a perturbation term $p^{\prime}$; $\rho$ is the seawater density, $\rho_{c}$ is the reference density (1035 kg m^-3 in our simulations), and $\rho^{\prime}$ is the density variation relative to the reference value; $g$ is the gravitational acceleration; $\lambda$ is the longitude; $\varphi$ is the latitude; $z$ is the vertical distance from the surface and is negative for seawater below the surface; $a=R+z$ is the distance from the center of the Earth, where $R$ is the planetary radius; generally $z\ll R$ and the approximation $a\approx R$ is used in models; ${D\over Dt}={\partial\over{\partial t}}+\vec{u}\cdot{\nabla}_{h}+w\,{\partial\over\partial z}$ is the total derivative, where ${\nabla}_{h}={1\over{acos\varphi}}{\partial\,\over\partial\lambda}\,\vec{i}+{1\over a}{\partial\,\over{\partial\varphi}}\,\vec{j}$ ( $\vec{i}$ and $\vec{j}$ are the zonal and meridional unit vectors, respectively); $A_{h}$ and $A_{z}$ are horizontal and vertical viscosities; $k_{h}$ and $k_{v}$ are horizontal and vertical diffusivities for potential temperature; $k_{hS}$ and $k_{vS}$ are horizontal and vertical diffusivities for salinity; the Coriolis parameter ($f$) is defined as $2\Omega sin\varphi$, where $\Omega$ represents the planetary rotation rate; $F_{u}$, $F_{v}$, $F_{\theta}$, and $F_{S}$ are zonal wind stress, meridional wind stress, sea surface temperature forcing, and surface salinity forcing, respectively.

For simulations with shallow oceans, $g$ is set to be constant throughout the entire ocean. While for the simulations with a 40-km-deep ocean, the depth dependence of $g$ may be needed to be considered (Figure 1(a)), given by:

where $G$ is the gravitational constant, $M$ is the planetary mass, and $g_{0}$ is the gravitational acceleration at the sea surface and is set to be 9.81 m s^-2. It should be noted that gravity increases with depth when the core is significantly denser than water such as for terrestrial planets, and gravity should decrease with depth if the planet is homogeneous in density (Dragoni, 2020).

In the control experiment, most of the planetary and oceanic parameters resemble present-day Earth (Table 1). By default, $R$ is 6371 km, $\Omega$ is 7.27 $\times$ 10^-5 s^-1, $g_{0}$ is 9.81 m s^-2, and the ocean depth is 40 km. $A_{h}$ and $A_{z}$ are 1.2 $\times$ 10⁶ and 10^-3 m² s^-1, respectively. $k_{v}$ is 10^-5 m² s^-1 in the interior of the ocean and 10^-4 m² s^-1 at the sea surface due to the effect of wind stress mixing. $k_{h}$ is zero and the mixing of tracer properties along isopycnal is represented by the Gent-McWilliams scheme (GM, Gent & Mcwilliams 1990) with the Redi eddy parameterization (Redi, 1982). In this study, the surface salinity forcing is not considered. Thus, for simplicity $S$ is set to be 34.7 g kg^-1 everywhere and will not evolve with time, and both $k_{hS}$ and $k_{vS}$ are not used.

The equation of state (EoS) reflects the dependence of density on potential temperature, salinity, and pressure. In MITgcm, the nonlinear EoS ‘JMD95Z’ (Jackett & Mcdougall, 1995) is relatively accurate when the ocean depth is shallower than $\sim$10 km. Below 10 km, the accuracy of ‘JMD95Z’ is less. To calculate the density in a 40-km-deep ocean, ‘JMD95Z’ is modified following the ‘IAPWS-95’ formulation (Wagner & Pruß, 2002) that covers a validity range for temperatures from 0 ^∘C to 1000 ^∘C and pressures up to 10⁹ Pa. After modification, the accuracy of density is better guaranteed, as shown in Figure 1(b). The detailed relationship between the density of seawater and salinity, potential temperature, and pressure under the nonlinear EoS, ‘JMD95Z’, and modified ‘JMD95Z’, is given in Appendix A.

The ocean circulation is forced by the atmosphere through the wind stresses and sea surface temperatures, but the full interaction between the atmosphere and ocean is not considered in this study. Zonal and meridional wind forcings, $F_{u}$ and $F_{v}$, are given by $\tau_{x}\over{\rho_{c}\Delta z_{s}}$ and ${\tau_{y}\over{\rho_{c}\Delta z_{s}}}$, respectively, where $\tau_{x}$ and $\tau_{y}$ are zonally-symmetric zonal and meridional wind stresses, and $\Delta z_{s}$ is the depth of the surface layer of the ocean. By default, the wind stress resembling the present-day Earth is imposed (solid lines in Figure 2(a) & (b)). The sea surface temperatures are set to be restored to a zonally-symmetric profile ($\theta_{\ast}$) similar to current Earth (Figure 2(c)) with a relaxation timescale ($\tau_{\theta}$) being equal to 30 Earth days, given as $F_{\theta}=-{{1\over\tau_{\theta}}(\theta-\theta^{\ast})}$. Both the wind stress and sea surface temperature forcings are applied only in the surface layer of the model.

We perform series of sensitivity simulations to investigate the depth of the thermocline under different planetary and oceanic parameters, including ocean depth, EoS for seawater, depth dependence of gravitational acceleration, parameterization scheme for mesoscale eddies, viscosity, wind stress, vertical diffusivity, and planetary rotation rate. We briefly outline our procedures for these simulations below and Table 2 summarizes the corresponding setup and the range for each parameter.

By default, the ocean depth is set to be 40 km. In order to compare that with shallow oceans, we do two experiments with ocean depths of 5 km and 10 km. In the shallow ocean simulations, the variation of gravitational acceleration with depth is not considered.

Wind forcing and vertical mixing are two main factors that make the ocean circulation driven by surface density/heat forcing not be restricted within a thin layer at the surface but reach deeper (Vallis, 2019). The magnitude of wind stresses ranging from 0.5 to 4 times that in the control experiment is investigated. Generally, the magnitude of wind stresses is determined by atmospheric density, surface wind speed, and drag coefficient (Stewart, 2008). Olson et al. (2020) showed that due to the opposite influences of surface pressure on atmospheric density and on surface wind speed, the wind stresses only increase by a factor of two when the surface pressure is increased by a factor of ten. Also, a case with stronger and wider easterly winds and equatorward winds in the tropics is carried out (see the dashed lines in Figure 2(a) & (b)), which considers wider Hadley cells induced possibly by a slower planetary rotation rate.

The strength of mixing is quantified as vertical eddy diffusivity, which can range from about 10^-5 (background level in the ocean interior) to 10^-3 m² s^-1 (enhanced mixing due to such as bottom topography) on Earth (Waterhouse et al., 2014). Si et al. (2021) use simple scalings and obtain the possible strength of vertical mixing on asynchronous rotating habitable exoplanets around M-dwarfs, which is roughly 100 times the mean value on Earth. Interior diffusivities varying between 10^-6 and 2 $\times$ 10^-3 m² s^-1 are tested in this study. When the mixing is strong, the equilibrated system is no longer wind-dominated but mixing-dominated.

Planetary rotation rates of exoplanets may range from several hours to several hundred Earth days (Akeson et al., 2013; Impey, 2013). The rotation rate can directly affect the strength of Ekman pumping and subduction and then affect the wind-driven circulation (Vallis, 2019). Planetary rotation rates ranging from 1/30 to 24 times that of Earth under both low-diffusivity (10^-5 m² s^-1) and high-diffusivity (10^-3 m² s^-1) cases are examined in this study.

The GM scheme is always used to parameterize the effect of meso-scale eddies on tracer transports when using a non eddy-resolving ocean circulation model. Studies on Earth show that the GM scheme can effectively resolve the effect of meso-scale eddies and exert an evident impact on the temperature and the thermocline (Danabasoglu & McWilliams, 1995). In mid-latitudes, where the meridional temperature gradients are large and baroclinic eddy activities are strong, the usage of GM scheme instead of a specified horizontal diffusion can greatly reduce the temperature bias between models and observations. One simulation with GM scheme turned off is carried out.

To verify the influence of the modified EoS on the thermocline depth, we also run two simulations using a linear EoS with a thermal expansion coefficient of 2$\times$10^-4 m² s^-1 or using the original ‘JMD95Z’ (Figure 1(b)). The effect of excluding the depth dependence of gravitational acceleration in 40-km-deep ocean is also tested. The horizontal and vertical viscosities decreased by a factor of 2 and 10 at the same time are also tested due to their possible effects on the thermocline depth.

Due to the zonal symmetry of the large-scale circulation on water-rich exoplanets, simulation results using zonally-symmetric two-dimensional (2D, y–z) models are almost the same to that of three-dimensional (3D) models (see Section 3 below). Thus, 2D models are used in most of our simulations due to computational resource limitations.

We do not consider any topography at the bottom of the ocean, and a flat-bottomed ocean is used. But, a linear bottom drag is included at the sea floor with a linear coefficient of 10^-3 m² s^-1. Horizontal resolution of the model is 2.25${}^{\circ}\times 2.25^{\circ}$ in longitude and latitude, respectively. In the vertical, we use 70 unequally spaced levels for the 40-km-deep ocean. The vertical resolution increases nonlinearly from 20 m at the sea surface to 665 m at the sea bottom. When simulating the 5-km-deep ocean and 10-km-deep ocean, we employ 15 and 24 layers in the vertical, respectively. Our simulations are initialized from a state of rest and are integrated until a statistical equilibrium is reached. For most of our cases, the system achieves quasi-equilibrium within about 15,000 Earth years. If not mentioned, the results shown below are mean states obtained by averaging over the last 1000 years of each integration.

In this paper, we investigate the thermocline depth on water-rich exoplanets where the ocean depth can reach tens of to hundreds of kilometers and there is no any exposed continent, using the global ocean model MITgcm. Our main conclusions are as follows:

1). The depths of the thermocline in our simulations are restricted within the upper ocean in most cases. In all our simulations, the maximum depth that the thermocline can reach is about 10 km. Due to the lack of exposed continent, the scaling theories for the thermocline depths in both wind-dominated system and mixing-dominated system are a little different from that used on Earth. And, the scaling theories are roughly consistent with our numerical results.

2). The wind-influenced thermocline depth is mainly determined by the wind forcing and becomes deeper with larger wind stress curl. The influence of planetary rotation rate on the global-mean depth of wind-influenced thermocline is limited. Varying the equation of state and not considering the depth dependence of gravity acceleration have limited influences on the thermocline depth. The influence of turning off the GM scheme for mesoscale eddies is relatively large, which makes the thermocline reach a depth of $\sim$8 km from a depth of $\sim$2 km. The thermocline can become deeper as the horizontal and vertical viscosities decreases, but the change is small.

3). The mixing-influenced thermocline becomes deeper as the vertical eddy diffusivity or planetary rotation rate increases, but overall it can not reach the bottom of the ocean under proper parameters.

Due to the lack of exposed continents and coastal upwelling motions on water-rich exoplanets, the depth of the ocean overturning circulation plays an important role in transporting the nutrients and other materials like carbon upward from the deep ocean to the upper ocean, which then affects planetary habitability. In our simulations, there are non-zero vertical velocities throughout the entire ocean under different magnitudes of planetary rotation rate, wind stress, and diffusivity (Figure 6 and upper panels of Figure 12), which suggests that the ocean overturning circulation might reach the bottom of the ocean. The magnitudes of the mean vertical velocity in our simulations are comparable to or weaker than the magnitudes of the mean vertical velocity in the Earth’s oceans (Figure 12). Thus, the ocean overturning circulation on water-rich exoplanets might reach the bottom of the ocean and help transport nutrients back to the thermocline and to the surface, which is beneficial to the biological activity and habitability on water-rich exoplanets.

In this study, the effect of varying wind stresses is tested without considering the full interaction between the atmosphere and the ocean and the imposed wind stresses are steady, so ocean-atmosphere coupling model can be utilized to investigate the thermocline depth in future work. The highest vertical diffusivity for the interior ocean we tested is 10^-3 m² s^-1 or 100 times that on Earth, which is the estimated strength of vertical mixing on potentially habitable asynchronous rotating exoplanets around M-dwarfs (Si et al., 2021). However, this estimation is obtained based on simple scalings, and more rigorous calculations of the background vertical diffusivity using high-resolution ocean models are required for further investigations.

The possible local volcanic activity at the bottom of the ocean is not included in our simulations. On Earth, the geothermal heat added at the sea floor by local volcanic activity is an important source for regional mixing in the deep ocean (Huang, 1999), which could affect the overlying density fields and the stratification of the ocean and then influence the strength and depth of the meridional overturning circulation (Adcroft et al., 2001; Mullarney et al., 2006). Thus, the possible local volcanic activity might also exert an influence on the thermocline depth on water-rich exoplanets. What’s more, the existence of high-pressure ice at the sea floor, which is an important characteristic of water-rich exoplanets, is also not included in our study. The high-pressure ice can reduce the bottom friction, and the decreased friction could result in a decreased strength of vertical mixing, which is unfavorable for the thermocline to reach deeper. The existence of the high-pressure ice can also directly inhibits the nutrient exchange between the ocean and the solid surface. Detailed simulations considering the effects of local volcanic activity or high-pressure ice are left for further investigations.

For simplicity, the effect of salinity is not included in this study. Both salinity forcing and temperature forcing at the surface can induce density differences and then drive ocean overturning circulation (Vallis, 2019). And, it is the strength of vertical mixing that largely determines the thermocline depths under both thermal-driven circulation and salinity-driven circulation. In this aspect, the effect of temperature forcing is basically similar to the effect of salinity forcing, and the results of simulations here that only consider the density changes induced by temperature forcing could be used to salinity-forcing experiments.

In this study, our model has a relatively coarse horizontal resolution of 2.25 ^∘ $\times$ 2.25 ^∘ and can not resolve mesoscale or smaller eddies explicitly. A high-resolution or eddy-resolving ocean model is required to confirm the reliability of our results. Ashkenazy (2017) used MITgcm with a high horizontal resolution of 100 m in a 10 km $\times$ 10 km domain with an ocean depth of 1.2 km to investigate the energy transfer from the surface ocean to the deep ocean. His results showed that the wind-induced currents are restricted within the top shallow layers except when the wind stress is temporally periodic with the Coriolis frequency, i.e., the wind forcing is resonating with the Coriolis force. We reproduce his two experiments but employ a 10-km-deep ocean and obtain the same result: the wind-induced kinetic energy is restricted within the surface layer under steady wind forcing but can be transferred to the deep ocean via the resonance between the wind forcing and the Coriolis force (figures not shown). However, this resonance occurs only in very limited region(s) where the frequency of the wind stress is exactly equal to the frequency of the Coriolis force, so it can not have a significant effect on the conclusions shown in this study.