Thermocline Depth on Water-rich Exoplanets

Water-rich exoplanets (alternatively called ocean planets) are 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 to hundreds of kilometers, which is much deeper than the global-mean depth of the ocean on Earth, ∼3 km (Stewart 2008). Due to the deep ocean, the pressure at the bottom of the ocean can be extremely high, which results in the occurrence of high-pressure ice at the seafloor. The existence of water-rich exoplanets and their possible bulk compositions has not yet been demonstrated directly by observations but is 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 better at resisting changes in external forcings (Hayworth & Foley 2020). Krissansen-Totton et al. (2021) suggested that the high pressure at the seafloor 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 increasing temperature (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 are no subtropical or subpolar gyres, boundary currents, and coastal upwelling motions on water-rich exoplanets owing 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 in the layers above or below (Pedlosky 2006; Stewart 2008; Fiedler 2010). The thermocline on Earth can be divided into two kinds: one is the main thermocline, and the other is the seasonal thermocline that varies with the seasons (Stewart 2008). The thermocline separates the warm mixed layer (turbulent and with 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 interregion mixing, oxygen content below the thermocline rapidly depletes with increasing depth, 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, the vertical density gradient within the thermocline is also larger than the upper mixed layer and the lower layer.

Here we first explore the depth of the thermocline under different ocean depths. Then, the sensitivity of the thermocline depth in wind-dominated systems and mixing-dominated systems 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, 1997b). 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 the 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

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{Du}}{{Dt}}-\displaystyle \frac{{uv}}{a}\tan \varphi -{fv}+\displaystyle \frac{1}{{\rho }_{c}a\cos \varphi }\displaystyle \frac{\partial {p}^{{\prime} }}{\partial \lambda }\\ -\,{{\rm{\nabla }}}_{h}\cdot ({A}_{h}{{\rm{\nabla }}}_{h}u)-\displaystyle \frac{\partial }{\partial z}\left({A}_{z}\displaystyle \frac{\partial u}{\partial z}\right)={F}_{u},\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{Dv}}{{Dt}}+\displaystyle \frac{{u}^{2}}{a}\tan \varphi +{fu}+\displaystyle \frac{1}{{\rho }_{c}a}\displaystyle \frac{\partial {p}^{{\prime} }}{\partial \varphi }\\ -\,{{\rm{\nabla }}}_{h}\cdot ({A}_{h}{{\rm{\nabla }}}_{h}v)-\displaystyle \frac{\partial }{\partial z}\left({A}_{z}\displaystyle \frac{\partial v}{\partial z}\right)={F}_{v},\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {p}^{{\prime} }}{\partial z}=-\rho ^{\prime} g,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{h}\cdot {\boldsymbol{u}}+\displaystyle \frac{\partial w}{\partial z}=0,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\rho ^{\prime} =\rho (\theta ,S,{p}_{0}(z))-{\rho }_{c},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\displaystyle \frac{D\theta }{{Dt}}-{{\rm{\nabla }}}_{h}\cdot ({k}_{h}{{\rm{\nabla }}}_{h}\theta )-\displaystyle \frac{\partial }{\partial z}\left({k}_{v}\displaystyle \frac{\partial \theta }{\partial z}\right)={F}_{\theta },\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{DS}}{{Dt}}-{{\rm{\nabla }}}_{h}\cdot ({k}_{{hS}}{{\rm{\nabla }}}_{h}S)-\displaystyle \frac{\partial }{\partial z}\left({k}_{{vS}}\displaystyle \frac{\partial S}{\partial z}\right)={F}_{S},\end{eqnarray} \tag{ 7 }$$

where u is the horizontal flow vector and u and v are its zonal and meridional components, respectively; w is the vertical current speed; θ and S are the potential temperature and salinity, respectively; the pressure field is splitted into a reference function of height p₀(z) and a perturbation term p'; ρ is the seawater density, ρ_c is the reference density (1035 kg m⁻³ in our simulations), and ρ' is the density variation relative to the reference value; g is the gravitational acceleration; λ is the longitude; φ 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 Earth, where R is the planetary radius; generally z ≪ R and the approximation a ≈ R is used in models; $\displaystyle \frac{D}{{Dt}}=\displaystyle \frac{\partial }{\partial t}+{\boldsymbol{u}}\cdot {{\rm{\nabla }}}_{h}+w\,\displaystyle \frac{\partial }{\partial z}$ is the total derivative, where ${{\rm{\nabla }}}_{h}=\displaystyle \frac{1}{a\cos \varphi }\displaystyle \frac{\partial \,}{\partial \lambda }\,{\boldsymbol{i}}+\displaystyle \frac{1}{a}\displaystyle \frac{\partial \,}{\partial \varphi }\,{\boldsymbol{j}}$ ( i and 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Ωsin φ, where Ω represents the planetary rotation rate; and F_u , F_v , F_θ , 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

$$\begin{eqnarray}&&g(z)=\displaystyle \frac{{GM}}{{\left(R+z\right)}^{2}}={g}_{0}{\left(\displaystyle \frac{R}{R+z}\right)}^{2},\end{eqnarray} \tag{ 8 }$$

where G is the gravitational constant, M is the planetary mass, and g₀ is the gravitational acceleration at the sea surface and is set to be 9.81 m s⁻². 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).

Zoom In Zoom Out Reset image size

In the control experiment, most of the planetary and oceanic parameters resemble present-day Earth (Table 1). By default, R is 6371 km, Ω is 7.27 × 10⁻⁵ s⁻¹, g₀ is 9.81 m s⁻², and the ocean depth is 40 km. A_h and A_z are 1.2 × 10⁶ m² s⁻¹ and 10⁻³ m² s⁻¹, respectively. Parameter k_v is 10⁻⁵ m² s⁻¹ in the interior of the ocean and 10⁻⁴ m² s⁻¹ at the sea surface owing to the effect of wind stress mixing. Parameter k_h is zero, and the mixing of tracer properties along isopycnal is represented by the Gent−McWilliams (GM) scheme (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⁻¹ everywhere and will not evolve with time, and neither k_hS nor k_vS is 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 ∼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 $\displaystyle \frac{{{\tau }}_{x}}{{{\rho }}_{c}{\rm{\Delta }}{z}_{s}}$ and $\displaystyle \frac{{\tau }_{y}}{{\rho }_{c}{\rm{\Delta }}{z}_{s}}$, respectively, where τ_x and τ_y are zonally symmetric zonal and meridional wind stresses, respectively, and Δ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 Figures 2(a) and (b)). The sea surface temperatures are set to be restored to a zonally symmetric profile (θ_*) similar to current Earth (Figure 2(c)), with a relaxation timescale (τ_θ ) equal to 30 Earth days, given as ${F}_{\theta }=-\displaystyle \frac{1}{{\tau }_{\theta }}(\theta -\theta * )$. Both the wind stress and sea surface temperature forcings are applied only in the surface layer of the model.

Zoom In Zoom Out Reset image size

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 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 by which the ocean circulation driven by surface density/heat forcing is not restricted within a thin layer at the surface but reaches 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 2 when the surface pressure is increased by a factor of 10. Also, a case with stronger and wider easterly winds and equatorward winds in the tropics is carried out (see the dashed lines in Figures 2(a) and (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⁻⁵ (background level in the ocean interior) to 10⁻³ m² s⁻¹ (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⁻⁶ m² s⁻¹ and 2 × 10⁻³ m² s⁻¹ 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⁻⁵ m² s⁻¹) and high-diffusivity (10⁻³ m² s⁻¹) cases are examined in this study.

The GM scheme is always used to parameterize the effect of mesoscale 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 mesoscale eddies and exert an evident impact on the temperature and the thermocline (Danabasoglu & McWilliams 1995). In midlatitudes, where the meridional temperature gradients are large and baroclinic eddy activities are strong, the use of the GM scheme instead of a specified horizontal diffusion can greatly reduce the temperature bias between models and observations. One simulation with the 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 × 10⁻⁴ m² s⁻¹ or using the original "JMD95Z" (Figure 1(b)). The effect of excluding the depth dependence of gravitational acceleration in a 40 km deep ocean is also tested. Decreasing the horizontal and vertical viscosities by a factor of 2 and 10 at the same time is also tested owing 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 as that of three-dimensional (3D) models (see Section 3). Thus, 2D models are used in most of our simulations owing to computational resource limitations.

We do not consider any topography at the bottom of the ocean, and a flat-bottomed ocean is used. However, a linear bottom drag is included at the seafloor with a linear coefficient of 10⁻³ m² s⁻¹. The horizontal resolution of the model is 225 × 225 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 yr of each integration.

3.1. Thermocline Depth in Wind-dominated System

We find that the thermocline depth on water-rich exoplanets with deep oceans is as shallow as that on planets with shallow oceans (Figures 3 and 4). In the three cases with different ocean depths (5, 10, and 40 km), cold water upwells near the equator and ±60° latitudes, and warm water at the surface sinks near ±30° latitudes, which results in lower temperatures in upwelling regions and higher temperatures in downwelling regions in the upper ocean (Figure 3(a)). Thus, the temperature field exhibits a bowl-shaped pattern. The potential density field is similar to the potential temperature field (Figure 3(b)) owing to the spatially uniform salinity used in this study. In the meantime, there are poleward flows in the tropics and equatorward flows in the midlatitudes in the surface layer (Figure 3(c)). Thus, the wind-driven meridional overturning circulation is clockwise in the tropics and anticlockwise in the midlatitudes. Ekman pumping and subduction induced by the imposed wind stresses can account for the equilibrated wind-driven overturning circulation (Vallis 2019), and the vertical velocity at the bottom of the Ekman layer can be given as

$$\begin{eqnarray}&&{W}_{{\rm{E}}}=\displaystyle \frac{1}{{\rho }_{c}}{\mathrm{curl}}_{z}\displaystyle \frac{\tau }{f},\end{eqnarray} \tag{ 9 }$$

where τ is the imposed wind stress and curl_z takes its curl in the vertical. Equation (9) shows that wind-driven Ekman pumping and subduction induce upwelling of cold water where the wind stress curl is positive in the northern hemisphere (negative in the southern hemisphere) and downwelling of warm water where the wind stress curl is negative in the northern hemisphere (positive in the southern hemisphere). The simulation results (Figures 3 and 5) are consistent with Equation (9).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The meridional currents at the sea surface are in the balance between the Coriolis force and the zonal wind stresses (Appendix B), so the meridional flow must be poleward when the zonal wind stress is easterly and equatorward when the zonal wind stress is westerly in the northern hemisphere (Figure 3(c)).

Equation (9) also suggests that there is no relationship between the ocean depth and the strength of Ekman pumping (subduction), i.e., the strength of wind-driven overturning circulation may be the same under different ocean depths, which can be found in Figure 3. The thermocline depth also exhibits little dependence on the ocean depth (Figures 4(a)–(c)). In all cases, the thermocline depth is restricted within the upper 2 km. Below 2 km, both the potential temperature and potential density are nearly uniform. Quantitative calculations based on the e-folding depth of potential density (Bryan 1987) show that the thermocline depths in the three cases with different ocean depths are nearly the same (Figure 4(c)). In the meridional direction, the thermocline is shallow near the equator, where the cold water is brought up and the vertical temperature (or potential density) gradient is relatively large, and the thermocline is deep at the midlatitudes, where the warm water sinks and the vertical temperature gradient is relatively small.

There is also no obvious difference in the sea surface height (SSH; Figure 3(e)) and zonal current fields (Figure 3(d)) among the three cases with different ocean depths. The latitudinal distribution of SSH is dominated by the effect of seawater temperature, i.e., SSH is high where the temperature of the underlying ocean is high and is low where the temperature of the underlying ocean is low. As a result, there is an SSH minimum near the equator owing to the upwelling of cold water and an SSH maximum near 30° because of the downwelling of warm water. The meridional gradient of SSH determines the meridional gradient of seawater pressure and thereby the speed of zonal currents through geostrophic balance both at the ocean surface and in the ocean interior (Figure B1 in Appendix B). Thus, the ocean flows are westward in the tropics and eastward in the midlatitudes (Figure 3(d), colors). With depth increasing, the speed of the zonal currents decreases, following thermal wind balance (contour lines in Figure 3(d); Vallis 2019).

As the ocean becomes deeper, the zonal flows become closer to the thermal wind balance, as the calculated thermal wind currents do not coincide as well with the zonal flows in the 5 km deep-ocean case as that in the 40 km deep-ocean case (Figure 3(d)). This is probably because the influence of the bottom drag on the zonal flows becomes weaker as the ocean becomes deeper.

The thermocline depth with varying ocean depths is also tested using a 3D configuration. Zonally symmetric external forcings, shown as the solid lines in Figure 2, are employed in the 3D simulations. Due to the axisymmetric external forcings and to the lack of exposed continents, the results of the 3D simulations are quite similar to the 2D simulation results (see Figure 4). It is worth noting that there is baroclinic instability within approximately 30°S–30°N in the 3D simulations, while there is not in the 2D simulations. Despite the existence of baroclinic eddy activities, its effect on potential temperature is quite limited. For details, see Appendix C. Given the consistency between the results of 3D and 2D simulations, we hereafter choose to use a 2D configuration to investigate the relationship between the depth of ocean circulation and other parameters in order to save computational sources.

3.1.1. The Effect of Varying Wind Stresses

The wind-driven overturning circulation becomes stronger and the wind-influenced thermocline reaches deeper as the magnitude of the wind stress curl increases (Figures 5(a)–(c)). As Equation (9) suggests, larger wind stress curl could induce stronger Ekman pumping and then result in stronger overturning circulation. Global-mean potential density profiles under different strengths of wind stress indicate that stronger wind forcings correspond to a deeper thermocline (Figure 5(a)), despite that the changes of the thermocline depth are highly latitude dependent (Figure 5(b)). With stronger wind forcings, the upwelling motion becomes stronger and then brings colder water to the upper ocean, which results in larger vertical temperature gradients and a shallower thermocline; meanwhile, the downwelling motion in midlatitudes is stronger and pumps more warm water deeper to the interior ocean, which results in weaker vertical temperature gradients and a deeper thermocline (Figure 5(b)). Overall, the global-mean depth of the thermocline increases with wind forcing, which extends from about 0.6 to 1.0 km when the wind forcing is increased by a factor of 4 (Figure 5(c)). The influence of varying the spatial pattern of imposed wind stress on the thermocline depth is limited (green square in Figure 5(c)), possibly due to the limited changes of the wind stress curl (Figures 2(a) and (b)).

A simple scaling for the depth of the wind-influenced thermocline can be derived from thermal wind balance and zonally symmetric mass continuity under the Boussinesq approximation (Welander 1971; Bryan 1987; Vallis 2000),

$$\begin{eqnarray}&&f\displaystyle \frac{\partial u}{\partial z}=-\displaystyle \frac{\partial b}{\partial y},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial v}{\partial y}+\displaystyle \frac{\partial w}{\partial z}=0,\end{eqnarray} \tag{ 11 }$$

where $b=-\displaystyle \frac{\rho ^{\prime} }{{\rho }_{c}}g$ is the buoyancy. Here $b\,=\,b^{\prime} (y,z,t)\,+\overline{b(z)}$ is composed of two parts: $\overline{b}$ is a reference field and only depends on depth; b' is the perturbation term that varies in both time and space. If considering a linear EOS and excluding the salinity effect, we have ρ = ρ_c (1 − αδT) and b' = αδTg, where α is the thermal expansion coefficient and δT is the temperature anomaly to a reference state. Thus, the scaling relationships given by Equations (10) and (11) are, respectively,

$$\begin{eqnarray}&&\displaystyle \frac{{fV}}{D}\sim g\alpha \displaystyle \frac{{\rm{\Delta }}T}{L},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\displaystyle \frac{V}{L}\sim \displaystyle \frac{{W}_{{\rm{E}}}}{D},\end{eqnarray} \tag{ 13 }$$

where V is the horizontal velocity scale, Ekman pumping velocity W_E scales the vertical velocity under wind forcing (Equation (9)), ΔT is the characteristic meridional temperature difference across the thermocline, and D and L are vertical and horizontal distance scales, respectively. These scaling relationships yield the depth of the wind-influenced thermocline:

$$\begin{eqnarray}&&D\sim {\left(\displaystyle \frac{{W}_{{\rm{E}}}{{fL}}^{2}}{g\alpha {\rm{\Delta }}T}\right)}^{1/2}.\end{eqnarray} \tag{ 14 }$$

Equation (14) is commonly used to scale the wind-influenced thermocline depth on Earth. But it should be noted that Equation (14) is obtained by assuming that the zonal flow speed and meridional flow speed are in the same magnitude (Equations (12) and (13)), which is always true for Earth (Vallis 2019). However, this assumption might not be necessarily true in our simulations owing to the lack of exposed continents (e.g., Figures 3(c) and (d)). In our simulations, U and V follow the scaling as

$$\begin{eqnarray}&&{fV}\sim {A}_{h}\displaystyle \frac{U}{{L}^{2}},\end{eqnarray} \tag{ 15 }$$

where A_h is the horizontal viscosity coefficient. The U–V relationship is given by the balance between the horizontal dissipation and the Coriolis force in the interior ocean when the horizontal viscosity is relatively strong (Appendix B),

$$\begin{eqnarray}&&-{fv}={A}_{h}\displaystyle \frac{{\partial }^{2}u}{\partial {y}^{2}}.\end{eqnarray} \tag{ 16 }$$

Further, Equations (12), (13), and (15) yield

$$\begin{eqnarray}&&D\sim {{fL}}^{2}{\left(\displaystyle \frac{{W}_{{\rm{E}}}}{{A}_{h}g\alpha {\rm{\Delta }}T}\right)}^{1/2}.\end{eqnarray} \tag{ 17 }$$

Equations (17) and (9) suggest that the wind-influenced thermocline depth is proportional to the square root of the magnitude of wind stress curl, which is roughly consistent with our results (Figure 5(c)).

3.1.2. The Effect of Varying Planetary Rotation Rates

Equation (17) indicates that the wind-influenced thermocline depth increases with planetary rotation rate and scales as Ω^1/2, which is roughly consistent with the 30°S–30°N mean values in our simulations (red points and red line in Figure 5(f)). As planetary rotation rate increases, Ekman subduction near the equator becomes weaker (Equation (9) and Figure 6) and less cold water is brought up, which results in smaller vertical temperature gradients and a deeper thermocline in the low-latitude regions (Figure 5(e)). Thus, within the main upwelling regions of 30°S–30°N, the mean thermocline depth increases with planetary rotation rate, especially for slowly rotating cases.

Zoom In Zoom Out Reset image size

With increasing planetary rotation rate, both the upwelling motion near the equator and ±60° and the downwelling motion near ±30° become weaker. Thus, the thermocline depth becomes deeper in the upwelling regions but shallower in the downwelling regions under larger rotation rates (Figure 5(e)). At the same time, the latitudinal width of the upwelling zone at low latitudes decreases, and that of the downwelling regions extends equatorward at the same time (Figures 5(e) and 6). The narrowing of the upwelling regions and the widening of the downwelling regions under larger rotation rates suggest that the opposite change of thermocline depth may cancel out with each other and makes no net contribution to the global-mean value of the thermocline depth as shown in Figures 5(d) and 5(f). In addition, the mean thermocline depth of 30°S–30°N also increases at a slower rate under high rotation rate cases owing to the extended downwelling motion toward the equator (red points in Figure 5(f)).

3.1.3. Other Parameters

The sensitivity of the wind-influenced thermocline depth to other parameters is also examined (Figure 7). In the control case, the variation of gravitational acceleration with depth is considered and the EOS used is the modified "JMD95Z" (Table 2). Neglecting the depth dependence of gravitational acceleration (Figures 7(a) and (e)) or varying the EOS (Figures 7(b) and (f)) does not result in an obvious difference in the thermocline depth. This is probably because the wind-influenced thermocline is always confined within the upper ocean. However, turning off the GM scheme has a significant impact on the thermocline depth. As Figures 7(c) and (g) show, the ocean becomes warmer and the thermocline reaches deeper when the GM scheme is turned off, extending from a depth of ∼2 to ∼8 km. This result is consistent with previous studies, for example, Danabasoglu & McWilliams (1995), who suggested that the thermocline is more diffused than observed if a specified horizontal diffusion rather than the GM scheme is used. We also test the sensitivity of the wind-influenced thermocline depth to viscosity. We decrease both horizontal and vertical viscosity coefficients by a factor of 2 and 10, respectively, and find that the kinetic energy increases and the thermocline depth becomes deeper, but the variation is limited. For instance, the thermocline deepens from a depth of ∼0.5 to ∼2 km when the viscosity is decreased by a factor of 10 (Figures 7(d) and (h)).

Zoom In Zoom Out Reset image size

3.2. Thermocline Depth in Mixing-dominated System

Deep-ocean circulation (also called thermohaline circulation) is mainly driven by turbulent mixing (e.g., Vallis 2019), and the mixing-influenced thermocline depth on Earth is demonstrated to be sensitive to the vertical (or more precisely, diapycnal) eddy diffusivity (Bryan 1987; Hu 1996; Vallis 2000). To examine the relationship between vertical mixing and the thermocline depth on water-rich exoplanets, we design three series of simulations: upper-ocean enhanced mixing, vertically uniform mixing, and deep-ocean enhanced mixing. In the upper-ocean enhanced mixing cases, the vertical mixing with a diffusivity of 10⁻³ m² s⁻¹ is allowed in the upper 1 and 10 km, respectively (see blue and red lines in Figure 8(a)); for comparison, a case in which a strong, uniform vertical mixing throughout the whole ocean is also carried out (black line in Figure 8(a)). The effect of varying the strength of vertical mixing is also examined in the vertically uniform mixing scenarios, and the diffusivity ranging from 10⁻⁶ m² s^–1 to 2 × 10⁻³ m² s⁻¹ is tested. In the deep-ocean enhanced mixing cases, the vertical diffusivity is enhanced in the bottom 10 km of the ocean by two orders from the background value (Figure 8(b)). This series of cases is designed to mimic the effect of bottom topography on enhancing the mixing in the deep ocean (St. Laurent et al. 2002; Nikurashin & Ferrari 2013). Note that in the real oceans of Earth regions of enhanced mixing are approximately between the seafloor and above 0.5 or 1 km (St. Laurent et al. 2001; Waterhouse et al. 2014), whereas in our simulations we extend this decay scale to 10 km to more clearly see the effect of bottom topography if it exists. The effect of varying planetary rotation rates on the thermocline depth under high diffusivity is also examined in the vertically uniform mixing scenarios. The planetary rotation rate tested ranges from 1/30 to 8 times Earth's rotation rate. Again, for simplicity, salinity is set to be uniform everywhere, and only the density forcing induced by temperature difference is included in these simulations.

Zoom In Zoom Out Reset image size

There is no stratification where there is no mixing (Figure 9). Without wind forcing, the distribution of the equilibrated sea surface temperatures with latitude is similar to the imposed sea surface temperature forcing (see Figure 2(c)). In the ocean interior, the vertical mixing acts to diffuse heat from the upper, warmer layer downward to the lower, colder layer. Thus, the isotherms (isopycnals) extend downward to the ocean interior (Figures 9(a) and (b)). When the vertical mixing is restricted within the upper 1 or 10 km, the isotherms (isopycnals) are also forced to be restricted in the upper 1 or 10 km, respectively. Below those critical layers, the ocean is filled with the downwelling cold water from high latitudes and the potential density is vertically uniform (Figure 9(e)). When the strong vertical mixing is allowed throughout the entire ocean, the heat is diffused downward freely (right column of Figure 9), even though the isotherms extend to only about 10 km and cannot reach the bottom of the ocean (40 km) when the vertical diffusivity is as large as 10⁻³ m² s⁻¹, which is nearly the largest diffusivity found in small local regions of Earth's oceans (Waterhouse et al. 2014).

Zoom In Zoom Out Reset image size

The spatial distribution of SSH with latitude is also similar to that of the imposed sea surface temperature forcing (Figure 9(d)). However, there is no SSH minimum near the equator as was found in the wind-driven circulation cases (Figure 3(e)); this is due to the lack of Ekman subduction and upwelling of cold water there. The meridional gradient of SSH increases as the vertical mixing acts deeper in the ocean (Figure 9(d)). This is because the meridional gradients of sea surface temperature are also diffused deeper when the heat is diffused downward deeper from ocean surface to ocean interior. The zonal currents at the sea surface are in geostrophic balance and are directly determined by the meridional gradients of SSH. Thus, the surface zonal currents become stronger as the vertical mixing is allowed to act deeper (Figure 9(c)). The zonal currents at the sea surface are eastward in the upper ocean, and there are two maxima near ±20° latitudes owing to the largest SSH gradients there.

The thermocline depth becomes deeper with stronger mixing (Figures 10(a)–(c)). When the diffusivity is low, the thermocline is shallow everywhere. As the vertical diffusivity increases, heat from the upper warm water is diffused downward deeper to the interior ocean, which weakens the vertical stratification and then deepens the thermocline at all latitudes (Figures 10(a) and (b)). It is also worth noting that the thermocline depth exhibits more obvious variation with latitude under high-diffusivity cases (e.g., higher than 10⁻⁵ m² s⁻¹). In addition, the thermocline is deep in high latitudes owing to the cold water sinking and weak stratification there, and it is shallower in low latitudes. Moreover, the slope of the global-mean thermocline depth as a function of diffusivity is approximately 1/3 in the double logarithm coordinates (Figure 10(c)). This can be explained using simple scaling theories (Welander 1971; Vallis 2000). Considering that there is no wind forcing and the mixing is strong, we scale the depth of the mixing-influenced thermocline using the same thermal wind balance (Equation (10)), mass continuity (Equation (11)), and the U–V relationship on water-rich exoplanets provided by the balance between the horizontal dissipation and the Coriolis force in the ocean interior (Equation (16)). What is more, the thermodynamic equation considering the effect of diffusive terms needs to be included in high-diffusivity cases, given as

$$\begin{eqnarray}&&\displaystyle \frac{{Db}^{\prime} }{{Dt}}+{{wN}}^{2}={k}_{v}\displaystyle \frac{{\partial }^{2}b}{\partial {z}^{2}},\end{eqnarray} \tag{ 18 }$$

where ${N}^{2}=\displaystyle \frac{d\overline{b}}{{dz}}$ is the vertical gradient of the reference buoyancy and represents the vertical stratification of the system. The thermodynamic equation obeys an advective–diffusive balance, so the scaling relationships of Equations (10), (11), (15), and (18) become

$$\begin{eqnarray}&&\frac{{fV}}{\delta }\,\sim \,g\alpha \frac{{\rm{\Delta }}T}{L},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\frac{V}{L}\,\sim \,\frac{w}{\delta },\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{fV}\,\sim \,{A}_{h}\frac{U}{{L}^{2}},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\frac{w{\rm{\Delta }}T}{\delta }\,\sim \,\frac{{k}_{v}{\rm{\Delta }}T}{{\delta }^{2}},\end{eqnarray} \tag{ 22 }$$

where δ is the depth of the mixing-influenced thermocline and w is no longer the imposed wind-driven Ekman velocity but is internally determined. Assuming that the vertical temperature/buoyancy variations across the thermocline are comparable to its meridional variations, the mixing-influenced thermocline depth is

$$\begin{eqnarray}&&\delta \sim {\left(\displaystyle \frac{{k}_{v}{f}^{2}{L}^{4}}{{A}_{h}g\alpha {\rm{\Delta }}T}\right)}^{1/3}.\end{eqnarray} \tag{ 23 }$$

This relationship suggests a scaling of k_v ^1/3, which is quite consistent with our numerical results (Figure 10(c)). Even though it is obtained from high-diffusivity assumptions, this scaling seems still applicable to low-diffusivity cases (e.g., when the diffusivity is lower than 10⁻⁵ m² s⁻¹).

Zoom In Zoom Out Reset image size

The thermocline depth is also sensitive to the planetary rotation rate under high diffusivity. As Equation (23) suggests, the thermocline becomes deeper as the planetary rotation rate increases and scales as Ω^2/3 under high-diffusivity conditions, which is roughly consistent with our global-mean results (Figures 10(d)–(f)). For example, the global-mean thermocline depth increases from ∼1 to ∼5 km when the planetary rotation rate increases from 1/10 Ω_e to 10 Ω_e , where Ω_e is Earth's rotation rate. The planetary rotation rate influences the mixing-influenced thermocline by affecting the internally determined upwelling velocity. Combining Equations (22) and (23), the vertical velocity is

$$\begin{eqnarray}&&w\sim {{k}_{v}}^{2/3}{\left(\displaystyle \frac{{A}_{h}g\alpha {\rm{\Delta }}T}{{f}^{2}{L}^{4}}\right)}^{1/3},\end{eqnarray} \tag{ 24 }$$

which becomes larger as the planetary rotation rate decreases. Larger upwelling velocity in the ocean interior under smaller rotation rate acts oppositely to the effect of diffusion and can result in larger vertical temperature gradients and a shallower thermocline, which is roughly consistent with the numerical results as shown in Figure 10(f).

In the third scenario, the bottom-enhanced mixing has a limited effect on the thermocline depth (Figure 11). Wind stress forcing is included in this group of experiments. The effect of wind-driven Ekman subduction and the consequent cold water upwelling in the tropics is still obvious in the temperature field when the vertical diffusivity is 10⁻⁵ m² s⁻¹ in the upper 30 km (Figure 11(a)). When the vertical diffusivity in the upper ocean is increased to 10⁻³ m² s⁻¹, the strong heat diffusion dominates the temperature field and the effect of wind-induced cold water upwelling in the tropics is no longer evident (Figure 11(b)). In the meantime, heat is diffused downward deeper to the interior ocean as the vertical diffusivity increases, which results in a deeper thermocline. However, the thermocline is still restricted within the upper 10 km when the diffusivity in the bottom 10 km of the ocean is increased to as large as 10⁻¹ m² s⁻¹ (Figure 11(b)). Again, results of 3D simulations are consistent with the results of the zonally symmetric 2D simulations (Figures 11(c) and (d)).

Zoom In Zoom Out Reset image size

The seafloor bathymetry, which is not explicitly included in our simulations, could affect not only the strength of mixing but also the viscosity at the bottom of the ocean. In our simulations, the bottom topography-induced friction is parameterized with a linear bottom drag at the seafloor. In addition, simple tests suggest that the influence of varying the magnitudes of the drag coefficient on the depth of the thermocline is small (figures not shown).

In this paper, we investigate the thermocline depth on water-rich exoplanets where the ocean depth can reach tens to hundreds of kilometers and there is no 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 systems and mixing-dominated systems are a little different from that used on Earth. In addition, 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 EOS and not considering the depth dependence of gravity acceleration have limited influence 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 ∼8 km from a depth of ∼2 km. The thermocline can become deeper as the horizontal and vertical viscosities decrease, 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 cannot 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 nonzero vertical velocities throughout the entire ocean under different magnitudes of planetary rotation rate, wind stress, and diffusivity (Figure 6 and top 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 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.

Zoom In Zoom Out Reset image size

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 an 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⁻³ m² s⁻¹, 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 seafloor 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 is more, the existence of high-pressure ice at the seafloor, 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 inhibit 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). In addition, 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 in salinity-forcing experiments.

In this study, our model has a relatively coarse horizontal resolution of 225 × 225 and cannot 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 × 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 cannot have a significant effect on the conclusions shown in this study.

We are grateful to Yonggang Liu, Yongyun Hu, Xiaozhou Ruan, and Ru Chen for their helpful discussions. J.Y. is supported by the National Natural Science Foundation of China (NSFC) under grants 42161144011 and 42075046.

The nonlinear EOS "JMD95Z" (Jackett & Mcdougall 1995) as a function of salinity S, potential temperature θ, and pressure p is written as

$$\begin{eqnarray}&&\rho (S,\theta ,p)=\displaystyle \frac{\rho (S,\theta ,0)}{1-p/K(S,\theta ,p)}.\end{eqnarray} \tag{ A1 }$$

Note that p is in bars (10⁵ Pa) here. In addition, ρ (S, θ, 0) is a 15-term equation in terms of S and θ and is written as

$$\begin{eqnarray}&&\begin{array}{l}\rho (S,\theta ,0)=999.842594+6.793952\times {10}^{-2}\theta -9.095290\\ \times \,{10}^{-3}{\theta }^{2}+1.001685\times {10}^{-4}{\theta }^{3}\\ -\,1.120083\times {10}^{-6}{\theta }^{4}+6.536332\times {10}^{-9}{\theta }^{5}+8.24493\\ \times \,{10}^{-1}S-4.0899\times {10}^{-3}\theta S\\ +\,7.6438\times {10}^{-5}{\theta }^{2}S-8.2467\times {10}^{-7}{\theta }^{3}S+5.3875\\ \times \,{10}^{-9}{\theta }^{4}S-5.72466\times {10}^{-3}{S}^{\tfrac{3}{2}}\\ +\,1.0227\times {10}^{-4}\theta {S}^{\tfrac{3}{2}}-1.6546S\\ \times \,{10}^{-6}{\theta }^{2}{S}^{\tfrac{3}{2}}+4.8314\times {10}^{-4}{S}^{2}.\end{array}\end{eqnarray} \tag{ A2 }$$

K(S, θ, p) is the secant bulk modulus and is a 26-term equation in terms of S, θ, and p, which is written as

$$\begin{eqnarray}&&\begin{array}{l}K(S,\theta ,p)=1.965933\times {10}^{4}+1.444304\times {10}^{2}\theta -1.706103\\ \times \,{10}^{0}{\theta }^{2}+9.648704\times {10}^{-3}{\theta }^{3}\\ -\,4.190253\times {10}^{-5}{\theta }^{4}+5.284855\times {10}^{1}S-3.101089\\ \times \,{10}^{-1}S\theta +6.283263\times {10}^{-3}S{\theta }^{2}\\ -\,5.084188\times {10}^{-5}S{\theta }^{3}+3.886640\times {10}^{-1}{S}^{\tfrac{3}{2}}+9.085835\\ \times \,{10}^{-3}{S}^{\tfrac{3}{2}}\theta -4.619924\times {10}^{-4}{S}^{\tfrac{3}{2}}{\theta }^{2}\\ +\,3.186519\times {10}^{0}p+2.212276\times {10}^{-2}p\theta -2.984642\\ \times \,{10}^{-4}p{\theta }^{2}+1.956415\times {10}^{-6}p{\theta }^{3}\\ +\,6.704388\times {10}^{-3}{pS}-1.847318\times {10}^{-4}{pS}\theta +2.059331\\ \times \,{10}^{-7}{pS}{\theta }^{2}+1.480266\times {10}^{-4}{{pS}}^{\tfrac{3}{2}}\\ +\,2.102898\times {10}^{-4}{p}^{2}-1.202016\times {10}^{-5}{p}^{2}\theta +1.394680\\ \times \,{10}^{-7}{p}^{2}{\theta }^{2}-2.040237\times {10}^{-6}{p}^{2}S\\ +\,6.128773\times {10}^{-8}{p}^{2}S\theta +6.207323\times {10}^{-10}{p}^{2}S{\theta }^{2}.\end{array}\end{eqnarray} \tag{ A3 }$$

In modified "JMD95Z," the coefficient in the K(S, θ, p) for the term p is adjusted from 3.186519 to 2.5608 following "IAPWS-95" data (Wagner & Pruß 2002) to present the pressure dependence of density of seawater in deep oceans more accurately.

The momentum budget at steady states of the zonnaly symmetric 2D simulations is shown here. Take the control case in Table 2 as an example. Due to the zonal symmetry of the model, the zonal pressure gradient force is zero. In the zonal direction (Equation (1)), the zonal wind stress and the Coriolis force induced by the meridional currents balance with each other at the sea surface (Figure B1(a)); the contribution from the dissipation term at the sea surface is significant only near the equator. However, it is the dissipation term that balances with the Coriolis force in the ocean interior (Figure B1(b)). Further diagnostics show that it is the horizontal dissipation that largely dominates the dissipation term and balances the Coriolis term, while the contribution from vertical dissipation is small (top panels of Figure B2). The magnitude of horizontal viscosity coefficient can significantly affect the zonal momentum budget. When the horizontal viscosity is decreased by three orders of magnitude, the contributions from the other two terms, advection+metric and vertical dissipation, become dominant (bottom panels of Figure B2). When the horizontal viscosity coefficient is low, the overturning circulation is determined by the combination of wind stress and vertical diffusion/dissipation. In the meridional direction (Equation (2)), both the ocean surface (Figure B1(c)) and ocean interior (Figure B1(d)) are in geostrophic balance, i.e., the meridional pressure gradient force balances the Coriolis force induced by zonal currents. The contribution of the meridional wind stress is limited. The zonal-mean momentum budget of 3D simulations is quite similar to that of zonally symmetric 2D models (figures not shown).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It can also be seen that the acceleration in the meridional direction is one or two orders of magnitude larger than the acceleration in the zonal direction. This is because the zonal current is nearly an order of magnitude larger than the meridional current, which determines the relative magnitudes of the Coriolis force.

According to the Charney–Stern–Pedlosky criterion, baroclinic instability might occur if the meridional gradient of potential vorticity Q_y (=$\beta -\tfrac{{\partial }^{2}u}{\partial {y}^{2}}$, where $\beta =\tfrac{\partial f}{\partial y}=\tfrac{2{\rm{\Omega }}\cos \varphi }{a}$) changes sign in the ocean interior or Q_y is in the opposite sign to the vertical shear of zonal current u_z at the upper boundary (Vallis 2019). In our 3D simulations, Q_y is dominated by β, which is positive everywhere and does not change sign (Figure C1(a)). However, u_z at the upper boundary is negative and in the opposite sign to Q_y within about −30° to 30° (Figure C1(b)). Thus, there should be baroclinic instability there, which can be seen from the zonal current anomaly field approximately within −30° to 30° both at the sea surface (Figure C1(c)) and in the interior ocean (figures not shown).

Zoom In Zoom Out Reset image size

With baroclinic eddy activities, compared to that of 2D simulations, the temperature in 3D simulations is slightly increased and the meridional temperature gradient is decreased near ±30° (the edge of the area where baroclinic instability occurs) in the surface layer of the model (bottom panels of Figure C1). However, this effect is limited, possibly due to the fixed temperature forcing at the surface. There is also an evident temperature decrease in the interior ocean near the equator due to stronger upwelling in the 3D simulation, which results in larger temperature gradients there.