Slab ocean model

General description

The slab ocean model is primarily based on the work of Codron (2012). This model uses the same horizontal grid as the GCM and is composed of two layers. The first layer (50 m depth) represents the surface mixed layer, where the exchanges with the atmosphere take place. The second layer (150 m depth) represents the deep ocean. Each oceanic grid point of the model is defined by the sea ice fraction & thickness and snow thickness.

In general for a slab ocean, the temperature at the ocean surface evolves as per the following equation: 
 $$\rho c_p H \frac{\partial T}{\partial t} = F_{surf} + Q_{flux}$$

The evolution of temperature is driven by surface fluxes ($$F_{surf} = F_{rad} + F_{sens} + F_{lat}$$) arising from the atmosphere model and a prescribed / parameterised flux ($$Q_{flux}$$) that can represent the effect of geothermal heating / ocean circulation (depends on the kind of planet that is being simulated). Both, the surface fluxes and the conductive fluxes coming from within interact with snow or ice that can be present at the ocean surface. Over each grid point, sea ice is represented by a uniform layer of depth $$H$$ and fractional area $$f$$. It may be covered by a layer of snow, which has its own heat capacity. The temperature at the base of the ice layer is always made to be equal to the freezing temperature of sea water $$T_o$$, for reasons shown below.

In its current form in the Generic-PCM, the slab ocean model handles two key processes, (a) Sea ice evolution, and (b) Heat transport by ocean circulation:

Sea ice evolution

Sea ice develops when the ocean temperature falls below the freezing point of water ($$T_0=–1.8°C$$, the freezing point of Earth's sea water; function of salinity) and melts when its temperature rises above freezing. The changes in ice extent and thickness are driven by energy conservation, such that the ocean temperature is kept at $$T_0$$ as long as ice is present. Specifically, if the temperature of the surface layer instantaneously falls below its freezing point, its temperature is set back to $$T_0$$, and the resulting energy difference is used to create ice. A framework of rules determines how these energy disparities can make ice thickness and coverage evolve. Conversely, if the ocean temperature becomes positive when ice is already present, its temperature is set back to $$T_0$$, and the resulting energy difference is used to melt ice. The same arguments also apply to the melting of snow if it has precipitated onto the ice.

Variation of bare sea ice albedo as a function of ice thickness for the VIS (250-690 nm; dark green line) and NIR (690-4000 nm; red line) spectral bands. The curves have been fit to long-term Antarctic sea ice observations made by Brandt et al. (2005)

Further, a planet's simulated climate (and consequently, its observability) is highly sensitive to the albedo parameterisations of sea ice and snow of the GCM. Charnay et al. (2013) addressed this by treating the albedo of bare sea ice as a function of ice thickness. The oceanic surface albedo ($$A$$) could be that of snow (if any) or bare sea ice. In the model, the bare sea ice albedo increases between the minimum and maximum albedo values as a function of ice thickness (see the light green curve in the figure on the right) and is given by:

$$A = A^\rm{max}_\rm{ice} – (A^\rm{max}_\rm{ice} – A^\rm{min}_\rm{ice} ) e^{\left ( \frac{–h_\rm{ice}}{h^0_\rm{ice}} \right )}$$

Here, we improve upon the work of Charnay et al. (2013) by instituting a spectral dependence of albedo, complementary to the existing ice thickness dependence. In the above equation, $$A^\rm{min}_\rm{ice}$$ is the minimal albedo (that of liquid water), equal to 0.06 irrespective of the incoming spectra. This is because the albedo of liquid water is almost constant between the Visible (VIS) and Near-Infrared (NIR) bands. $$A^\rm{max}_\rm{ice}$$ is the maximal albedo, which we computed by fitting the above equation to long-term albedo observations of Antarctic sea ice documented by Brandt et al. (2005; see the figure on the right) and summarised by Pedersen et al. (2009). With this, $$A^\rm{max}_\rm{ice}$$ is found to be 0.69 and 0.31 in the VIS (250-690 nm; dark green curve) and NIR (690-4000 nm; red curve) bands, respectively. Further, $$h_\rm{ice}$$ is the ice thickness (in m) and $$h^0_\rm{ice}$$ is the minimum threshold of ice, taken to be 0.3 m.

While the bare sea ice albedo is resolved up to 4000 nm, the new model only distinguishes between the VIS and NIR bands concerning radiation extinction with depth. This extinction assumes solar-type surface irradiance fractions in these bands, implying that radiation beyond a certain wavelength [(> 1200 nm ?!)] is not transmitted. This solar assumption is also seen in other exoplanet-based GCMs like ROCKE-3D and ExoCAM. It is expected to be revised later to incorporate sensitivity to different stellar spectral irradiances. Nevertheless, the new spectral dependency of albedo provides a reasonable initial approximation.

It should be noted that changes in the sea-ice model, or in parameters such as the snow albedo, can have considerably large impacts on the simulated climate. This happens through changes in the ice extent and the global temperature through the ice-albedo feedback. This effect is, however, largely independent of the heat transport schemes presented in this page.

Heat transport by ocean circulation

The transport of heat by the ocean circulation is given by four components:

1. The impact of subgrid-scale eddies is represented by horizontal diffusion, with a uniform diffusivity in both layers.
2. [SB WILL ADD FIGURE] The mean wind-driven circulation is computed by calculating the Ekman mass fluxes in the surface layer from the wind stress, and taking an opposite return flow at depth. We then use this mass flux to advect the slab temperature horizontally. If we have a divergence in the horizontal component of the transport, upwelling is induced (and therefore, temperature advection) between the two layers. Conversely, convergence implies downwelling. The integrated oceanic horizontal transport (mass flux) due to wind stress in the Ekman layer is given by $$ \mathbf{V_E} = - \frac{1}{f} \mathbf{k} \times \mathbf{\tau} $$, where $$\mathbf{\tau}$$ is the net surface wind stress and $$f$$ is the Coriolis parameter (orthogonal to the wind stress). This simplified model reproduces the global meridional oceanic heat transport quite closely compared to an Oceanic-GCM.
3. The two-layer nature of the ocean model also facilitates a convective adjustment. For instance, if the surface layer is colder than the deep layer (for e.g., in winter at mid-high latitudes), the vertical convective adjustment ensures that the heat stored at depth is then restored to the surface. In effect, this simulates denser colder water at the surface descending to the depths.
4. [SB WILL ADD FIGURE] We parameterise ocean eddies with the Gent-McWilliams scheme, which is based on oceanic isotherms. These act to mix temperatures horizontally. Figure XX represents a simplified cross-sectional view of the ocean horizontally (across latitude) and vertically (across depth). The top and bottom halves represent the mixed and deep layers of the ocean respectively. One example of an isotherm has been overplotted. The figure illustrates how the ocean surface temperature changes from warmer to cooler waters from the equator to the pole. Eddies act to reduce the slope of the isotherms, which results in high-latitude surface heating and low-latitude depth cooling. This in turn re-stratifies the vertical ocean layers as shown in the figure, which increases the temperature difference between the surface and deep ocean, making Ekman transport larger. The horizontal advective velocity is directly proportional to the slope of the isotherms.

Other ocean circulation features like density-driven / thermohaline circulation and horizontal gyres are not present in the model. While these features play an important role in shaping regional climates on modern Earth, their effect is weaker on the global average; gyres would be absent in the case of a global ocean.

The 2024 version

Compared to the old model, the new model has the following changes (non-exhaustive):

1. It can be used in parallel mode.
2. It has a more realistic description of sea ice creation and evolution - simultaneous surface, side and bottom melting / freezing depending on fluxes.
3. There is a spectral dependence of bare sea ice albedo as a function of ice thickness.
4. Snow has an effective heat capacity.
5. Snow has "weight"; it can sink an ice block if there is too much of it.
6. Snow can be blown off by wind.
7. The two-layer ocean facilitates convective adjustment.
8. Diffusion can follow the Gent-McWilliams scheme + Eddy diffusivity.
9. Since Ekman Transport doesn't work at the equator (Coriolis term = 0), we prescribe an alternative "Sverdrup scheme" for regions close to the equator which relies on direct friction-driven wind transport.

Variables of the slab ocean model which can be written as outputs in diagfi.nc:

1. tslab1, tslab2: temperature of the surface and deep ocean layers respectively (K)
2. pctsrf_sic: grid fraction of sea ice
3. sea_ice: mass of sea ice (kg/m$$^2$$)
4. tsea_ice: the temperature of the surface in contact with the atmosphere (sea ice or snow on top of sea ice) (K)
5. dt_hdiff1, dt_hdiff2: heating of the surface and deep ocean layers by horizontal diffusion (W/m$$^2$$)
6. dt_ekman1, dt_ekman2: heating of the surface and deep ocean layers by Ekman transport (W/m$$^2$$)
7. dt_gm1, dt_gm2: heating of the surface and deep ocean layers by Gent-McWilliams eddy transport (W/m$$^2$$)

Technical aspects

The slab ocean model is activated in callphys.def with:

ok_slab_ocean = .true.

There are a host of transport schemes that one can switch on/off in regard to oceanic heat transport. In general, if one wants to use all schemes of ocean heat transport, they can keep the settings detailed below. Note that the following can only be set to true if ok_slab_ocean = .true.

Ekman transport: (default=.false.)

slab_ekman = .true.

Ekman zonal advection: (default=.false. ; can only be true if slab_ekman = .true.)

slab_ekman_zonadv = .true.

Gent-McWilliams Scheme (can only be true if slab_ekman = .true.):

slab_gm = .true.

Horizontal diffusion (default=.false.):

slab_hdiff = .true.

Convective adjustment (0 - no, 1 - yes):

slab_cadj = 1

In the old implementation, tsea_ice was always the temperature of the ice surface. It didn’t matter if snow was on top of it or not because previously, snow had zero heat capacity. But now, since snow has a finite heat capacity, there can be conductive fluxes between the top layer of the sea ice and snow on top of it (if the latter exists). This necessitates two different temperature variables: (1), for the sea ice temperature if snow exists on top of it (tice, which is a global variable within ocean_slab_mod), and is arguably not useful in physiq_mod and therefore not defined in it), and (2), the temperature of the surface in contact with the atmosphere, which can either be sea ice or snow (tsea_ice, in physiq_mod, as it is essentially a proxy for surface temperature, which is useful for the physics to compute interactions between the ocean and atmosphere). To summarise, tice and tsea_ice are only identical if snow does not exist on top of the ice.

The timestep for updating the slab temperature and ice fraction is given by cpl_pas (in physics timesteps; default value = once every Earth day)

rnat is a standard diagfi output, which either represents 0 for oceanic grid points or 1 for continent grid points.

The model doesn't permit surface fluxes to enter the deep ocean. The latter is present in the model for exchanges and transport to take place.

The model also does not compute dynamical sea ice drift.

A word of caution: Ekman transport doesn’t work at the equator (Coriolis parameter = 0), we institute an alternative "Sverdrup scheme”, which relies on pressure-driven transport specifically at the equator. This scheme is also important for the case of slow-rotating planets when the Coriolis parameter is also low (?!).

QUESTION FOR FRANCIS: It seems that the Sverdrup scheme might only be beneficial for the special case of the Earth? Need to confirm.

The model follows energy conservation. However, if an upper limit on sea ice thickness is prescribed (for e.g., 10 m, for faster convergence), energy conservation can be violated. If one wants to purely conserve energy (at the cost of computation time), no limit on sea ice thickness can be prescribed. It is OK to prescribe an upper limit for the case of modern Earth (however, if your prescription is too small, it can lead to anomalous behaviour like stronger seasonal cycles). However, it might not be accurate to do this for more "exotic" cases like Snowball Earth and tidally locked ocean worlds.

Convergence - As far as the ocean is concerned, convergence timescales are governed by the thermal inertia and depth of the mixed layer. As an order of magnitude, convergence for Modern Earth with the dynamical slab ocean takes around 50 years.

---

SB: Latest edits (r3397, 17 Jul 2024) need to account for tice in newstart and start2archive also.