Difference between revisions of "Slab ocean model"

Revision as of 13:48, 2 April 2024

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. Further, the surface albedo ($$A$$) can be that of snow (if any) or that of bare ice 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, $$A$$ represents the albedo, $$A^\rm{max}_\rm{ice}$$ = 0.65 is the maximal albedo, $$A^\rm{min}_\rm{ice}$$ = 0.2 is the minimal albedo, $$h_\rm{ice}$$ is the ice thickness (in m) and $$h^0_\rm{ice}$$ = 0.5 m. The albedo over ice-free ocean is assumed to be 0.07. However, by definition, the Generic-PCM is meant to handle a more generalised case of planets. Since the albedo of a material (for e.g., ice or snow) is spectrally-dependent, the albedo of ice will be different around a star that is unlike the Sun. For instance, the ice/snow albedo climate feedback is significantly weaker for planets orbiting M-dwarf stars (for e.g., Proxima Centauri b) than for planets orbiting G-type stars like the Sun (see Joshi & Haberle [2012] and Shields [2014]). The Generic-PCM accounts for this by instituting a spectral dependence of albedo.

It should be noted that changes in the sea-ice model, or in parameters such as the snow albedo, can have 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.

Further, a planet's simulated climate and observability are highly sensitive to the albedo parameterisations of sea ice and snow in the GCM. We address this by treating the albedo of bare sea ice as a function of both, ice thickness and incoming spectra. Like in Charnay et al. 2013, the oceanic surface albedo ($$A$$) can be that of snow (if any) or that of bare ice. It increases between the minimum and maximum albedo values as a function of ice thickness 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. $$A^\rm{max}_\rm{ice}$$ is the maximal albedo, which can be 0.68 in the Visible (250-690 nm) regime and 0.31 in the NIR (690-4000 nm) regime. $$A^\rm{min}_\rm{ice}$$ is the minimal albedo (that of liquid water) and is equal to 0.06 irrespective of the incoming spectra, as the albedo of liquid water is almost constant between the VIS and NIR bands. We computed these albedo values by fitting the above equation to long-term albedo observations of Antarctic sea ice made by Brandt et al. 2005, in different spectral bands as a function

and is $$h_\rm{ice}$$ is the ice thickness (in m) and $$h^0_\rm{ice}$$ = 0.5 m.

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. Snow has an effective heat capacity.
4. Snow has "weight"; it can sink an ice block if there is too much of it.
5. Snow can be blown off by wind.
6. The two-layer ocean facilitates convective adjustment.
7. Diffusion can follow the Gent-McWilliams scheme + Eddy diffusivity.
8. 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: temperature of the sea ice surface (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

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.