GFDL - Geophysical Fluid Dynamics Laboratory

FMS Slab Ocean Model Technical Documentation

T. Knutson

Sept. 4, 2003

Revised Jan. 19, 2017

1. Introduction

In running climate experiments with FMS, the user has a choice of providing the atmospheric model with specified SST and sea ice boundary forcing or running in “coupled mode” with interactive ocean and sea ice models. For some purposes, the user may wish to couple the atmospheric model to a highly simplified “slab-type” ocean model rather than a full ocean general circulation model (GCM). In this document, the FMS slab ocean model (often also referred to as “mixed layer model” or “q-flux adjusted mixed layer model”) is described, along with some examples of its typical uses in climate sensitivity studies.

The remainder of this document is organized into the following sections:

  1. Background on Slab Model Experiments in Climate Studies

  2. Procedure for Performing a 2xCO2 Sensitivity Experiment

a. Restoring Run

b. Q-flux adjusted Control run

c. 2xCO2 perturbation run

  1. Description of the Model Equations

  2. Reconciliation of observed SST and sea ice climatologies

2. Background on Slab Model Experiments in Climate Studies

Running a climate experiment with a fully coupled model that includes a dynamical ocean GCM is a very computationally demanding task. Parts of the system, such as the deep ocean, require thousands of years of simulated evolution to equilibrate to a given radiative forcing. Also, before climate change experiments can be attempted, a long “control” integration with present-day or pre-industrial radiative forcing estimates is performed in order to demonstrate that the model is sufficiently stable, in terms of ocean circulation and thermal properties, to be useful for scientific experiments. In practice, full-coupled model climate change experiments are usually “transient” experiments in which the deep ocean does not fully equilibrate with the radiation perturbation.

On the other hand, if the equilibrium climate sensitivity is desired or if one wishes to preview the climate sensitivity of an atmosphere/land/sea ice model in a less computationally demanding setting, the typical approach is to use a slab model for the ocean component. In this case, the ocean is represented by a horizontal grid of points that are slabs of water of uniform specified depth (order 50 m) and salinity. These ocean points can have sea ice, which is modeled in FMS by the Sea Ice Simulator (SIS), documented in detail elsewhere. The FMS slab ocean code can be coupled either to a simple slab-type sea ice model without ice dynamics or to more physically complete sea ice representations that include dynamical effects such as ice advection by the winds. [As an aside, a future project will attempt to add currents to the slab ocean model for the purpose of advecting sea ice.]

The ocean temperature at each grid point in the slab ocean model is a predicted variable and is affected by heat exchange across the air-sea and ocean-sea ice interfaces. There is no direct communication between adjacent ocean grid points (i.e., no currents, temperature advection, diffusion, and so on), nor is there any explicit representation of the deep ocean.

The coupled system just described can be integrated forward in time and comes into equilibrium with a specified radiative forcing within a few decades. However, the many model simplifications–primarily the lack of ocean dynamics (advection), diffusion and convection, but also other model surface flux errors leads to a poor simulation of ocean temperature features such as the equatorial Pacific cold tongue that depend crucially on interactive ocean dynamics. Therefore in practice, a set of specified surface flux adjustments (commonly called “q-flux adjustments”) are developed and added to the slab model’s temperature tendency equation at each time step in order to maintain a seasonal cycle of ocean temperatures which is as close as possible to the present day conditions. The method for determining the q-flux adjustments is described in more detail in Section 3.

The coupled atmosphere/slab ocean/sea ice model with q-flux adjustments can then be intregrated for many decades with reasonable computational expense. Also for a climate change experiment, such as a doubling of atmospheric CO2, the perturbed model will equilibrate to the new forcing within a few decades, making it a feasible framework for performing benchmark experiments to document a given atmosphere/sea ice model’s “equilibrium” climate sensitivity. Owing to the need to run a fully coupled model (with full ocean GCM) for thousands of years to reach equilibrium, the “equilibrium sensitivity” of such fully coupled models has only rarely been computed.

3. Procedure for Performing a 2xCO2 Sensitivity Experiment

In this section, the procedure for performing a prototypical climate sensitivity dexperiment (i.e., a 2xCO2 equilibrium experiment) is described. The procedure consists of making three separate runs of the slab model: a restoring run, a q-flux adjusted control run (or 1xCO2 run), and a 2xCO2 perturbation run. The user must also select the type of sea ice model to couple to. If “SLAB_ICE” is set to “.true.” in the ice model namelist section, a simple slab-type ice model without ice dynamics is used. If “SLAB_ICE” is not set to “.true.” the default fully dynamical multi-partition sea ice model is used (see SIS documentation elsewhere in FMS for more details on this model.)

a. Restoring Run

A necessary step in performing a 2xCO2 sensitivity experiment is to prepare a control, or 1xCO2, experiment to serve as a baseline. However, before running the control experiment, however, the q-flux adjustments must be derived. These are derived in the “Restoring Run” of the slab mixed layer model, in which the SST and sea ice are restored toward climatological values. The ‘default’ climatology is part of FMS, and is currently based upon AMIP2 SST and sea ice climatologies. Restoring is turned on via namelist parameters for the SST and sea ice. It is possible to run the model restoring only SST or only sea ice, although this in not typically done. We typically use a 5-day restoring timescale for both SST and sea ice.Running in such a restoring mode, the slab model typically tracks the climatological seasonal cycle of SST and sea ice fairly well due to the strong restoring.

As the slab mixed layer model is integrated, the restoring heat fluxes are archived, and after a sufficient sampling period (we use 40 years, following an initial 2-year spin-up period which is discarded) a climatology of the restoring terms can be computed. This is obtained automatically when climatological versions of ocean_month and ice_month history files are generated. The q-flux adjustment is the sum of the sea-ice and SST restoring terms (all of which are in units of W/m2). The ice and SST restoring are combined because when the model is run in q-flux adjusted mode, the q-fluxes are applied to the mixed layer (ocean) and not directly to the sea ice, if present. The q-flux adjustment climatology can be computed, for example, in a Ferret session, since the ocean and ice climatology files are netcdf and can be read by Ferret, which then generates a netcdf file with the q-flux adjustment climatology.

b. Q-flux adjusted Control run

Once the q-flux adjustments have been derived, the slab model control run can begin. Restoring is turned off for SST and sea ice (separate namelist parameters), and q-flux adjustment is turned on (for the ocean model only). The q-flux adjustment climatology file, a netcdf file, must be supplied to the runscript. The model can then be integrated forward in time. The q-flux adjustments make up for the lack of ocean heat transport in the simple slab ocean and also correct for other errors or limitations in the model which lead to errors in the SST and sea ice simulations. However, the q-flux adjustments will not completely correct for such problems, so that some drift (order of 0.5 deg C globally for SST and some errors in simulated sea ice thickness and distribution) are to be expected. The user must make a subjective assessment of whether the simulation with the q-flux adjustments is giving a “good enough” simulation of the seasonal cycles of SST and sea ice or whether further work on the model is needed.

Examples of such further refinement to the model are discussed in some detail at this point. A very common problem with q-flux adjusted slab-type models is a runaway growth of sea ice thickness at some gridpoints in sea ice regions where the q-flux adjustments are negative (cooling) during at least part of the seasonal cycle. This effect can happen if the sea ice becomes unusually thick, due perhaps to climate variability in the model, and the ocean mixed layer beneath this ice becomes increasingly insulated from the overlying atmosphere. The negative q-flux adjustment will continue to work to take away heat, resulting in thicker sea ice, but the increasingly insulated mixed layer cannot receive heat as effectively from the atmosphere during seasons when it did so during the restoring runs. As a result, the mixed layer ocean becomes decoupled from the atmosphere and the sea ice continues to grow without bound.

To deal with this pathology of the q-flux adjusted slab ocean model, an artificial sea ice “lid” has been incorporated into the model. The sea ice lid can be turned on via a namelist parameter, and the maximum allowable sea ice thickness is also a namelist parameter specified by the user. We typically use 4 meters as the lid thickness, which is a compromise that allows most of the model gridpoints where the problem is occurring to equilibrate to the lid value in a reasonable (~40 years or less) time period. The sea ice lid operates as follows. If the sea ice thickness at a grid point exceeds the limit specified by the lid, the sea ice thickness is reduced to the lid value and the amount of heat needed to melt the amount of “shaved-off ice” is archived in the ice history variable QFLX_LIMIT_ICE. After the control model has reached a relatively stable equilibrium, including a relatively stable level of sea ice lid fluxes, the model is integrated for several decades to compute a stable baseline climatology, including a climatology of QFLX_LIMIT_ICE. In practice, we have typically integrated the control run for 100 years, allowing the first 50 years to be an equilibration period and the last 50 to be the relatively stable period for computing the control run climatology and the climatological sea ice lid fluxes.

A second example of a problem occurring at the control run stage that has, in our experience, required further refinements to the model is the “excessively cold eastern equatorial Pacific” problem. With that problem, the control run has drifted too cold by 8 degrees Celsius or more in the eastern equatorial Pacific. It was believed to be caused by a feedback loop involving excessive near-surface cloud, leading to cooler SSTs, which promoted even stronger low level cloud cover and so on. This problem was dealt with by imposing a lower limit “floor” on the vertical mixing coefficient at the lowest model level which had the effect of preventing the strongly stable situation in which the strong cloud cover feedback loop was being initiated. In later model versions, a new atmospheric boundary layer scheme was adopted based on the UK Met Office boundary layer; q-flux adjusted runs with this model to date have not exhibited the excessive cold tongue problem and thus have not needed any imposed lower limit on vertical mixing.

c. 2xCO2 perturbation run

Once the control run has been integrated long enough to derive a climatology of the sea ice lid fluxes (see previous section), the climate perturbation experiment can be started. Here we use a 2xCO2 experiment as an example. However, the same basic procedure applies for other types of perturbation experiments such as tropical deforestation or paleoclimate (e.g., last glacial maximum) experiments.

Before beginning the perturbation experiment, the q-flux adjustments are modified to include the heating contribution of the sea ice lid from the control run in order to conserve energy as the perturbation experiment must be run without an ice lid. This can be done, for example, in a Ferret session. The modified q-flux adjustment file must be identified in the runscript for the 2xCO2 run. The sea ice lid is turned off (namelist parameter) for the perturbation experiment, and the perturbation is specified (e.g., doubled CO2 concentration vs the control) in the namelist or script.

The perturbation model can then be integrated in time. We let the model equilibrate for about 50 years, keeping in mind that the slow parts of the slab model’s climate system, such as sea ice, take longer to equilibrate than land temperature, but not as long as the deep ocean for a fully coupled experiment (e.g., Stouffer et al., Clim. Dyn., 2003).

Once the full experiment is completed (typically 100 years of integration), the final 50 years are used to construct the perturbation run climatology, and the climate sensitivity, typically defined as the annual global mean difference in 2-m air temperature between the 2xCO2 and 1xCO2 experiments, can be computed.

4. Description of the Model Equations

In this section, the key model equations used in the mixed layer model are outlined and discussed. To facilitate comparison with the FMS model code, the notation and variable names used here is similar to that in the model subroutines.

The primary equation is the mixed layer temperature equation:

Ocean%t_surf = Ocean%t_surf + (dt_ocean * net_hflx/mlcp)


ocean%t_surf is the temperature of the mixed layer [K]
dt_ocean is the time step of the mixed layer model [sec]
net_hflx is the net heat flux into the mixed layer from the atmosphere and overlying sea ice [W/m2]
mlcp is the heat capacity of the mixed layer [J / m3 / deg C]
mlcp = mixed_layer_depth*4e6


mixed_layer_depth is the thickness of the mixed layer [meters, default of 50m]
4e6 is the assumed heat capacity of sea water in mixed layer [J / m2 / deg C]

The net heat flux is the sum of the heat flux heat flux into the mixed layer from the atmosphere and overlying sea ice plus any additional heat flux into the layer due to either restoring of the mixed layer temperature or a q-flux adjustment.

net_hflx(:,:) = hflx(:,:) + qflux_adj(:,:) + qflux_restore_sst(:,:)

The components of the heat flux from atmosphere or overlying sea ice are:

Hflx= Ice_boundary%sw_flux + Ice_boundary%lw_flux - (Ice_boundary%fprec + Ice_boundaryÊlving)*hlf - Ice_boundary%t_flux - Ice_boundary%q_flux*hlv

These components are the net downward flux of shortwave and longwave radiation, the heat flux due to the net flux of frozen precipitation directly from atmosphere (fprec) or land (calving), and the sensible and latent heat fluxes (t_flux and q_flux*hlv). With the latent heat of fusion (hlf) or vaporization (hlv) the units of these terms are all W/m2.

If SST restoring is turned on (a namelist parameter), it is computed by relaxing the temperature of the mixed layer toward the observed SST for that time of year and geographical location. The relaxation has a time scale defined by the namelist parameter sst_restore_timescale, which is input in units of days. The default value is 5 days for the SST.

qflux_restore_sst = (sst - Ocean%t_surf) * mlcp / (sst_restore_timescale * 86400)

In this equation, 86400 [sec/day] converts the restoring timescale from days to seconds.

If sea ice restoring is turned on (a namelist parameter) the model restores the sea ice enthalpy towards the observed sea ice enthalpy for that location and time of year.

For the full dynamical sea ice model, this means restoring towards the product of the observed thickness and the observed concentration. For the simpler slab ice, since the concentration is either zero or 100%, the restoring is towards the observed thickness.

This part of the mixed layer code is contained in the sea ice model routine ice_model.f90:

if (slab_ice) then
    heat_res_ice = -(LI*DI*Obs_h_ice(i,j)-sum(e2m)) &
    heat_res_ice = -(LI*DI*Obs_h_ice(i,j)*Obs_cn_ice(i,j,2)-sum(e2m)) &
end if


heat_res_ice is the restoring heat flux for the ice
Obs_h_ice is the observed sea ice thickness being restored to
Obs_cn_ice is the observed sea ice concentration being restored to (non-slab ice only)
Sum(e2m) is the enthalpy of the model sea ice
LI is the latent heat of fusion [J/kg]
DI is the density of ice [kg/m3]
dt_slow is the time step for the sea ice model
ice_restore_timescale is the relaxation time scale on which the model sea ice is restored, a namelist parameter with a default value of 50 days. 
86400 converts the restoring timescale from days to seconds

The enthalpy of the sea ice is computed for the simple “slab ice”-type sea ice model according to the expression given below. (A more complicated treatment for the non-slab ice, accounting for the sea ice concentration in various partitions of each ice grid cell, is contained in the model code, but not shown here.)

if (slab_ice) then
    e2m(2) = Ice%h_ice(i,j,2)*DI*LI


Ice%h_ice(i,j,2) is the thickness of the model's "slab" sea ice at each gridpoint.

Similar to ice restoring, if a limit is placed on the sea ice thickness to prevent unbounded growth, as discussed elsewhere in this document, the heat required to get rid of the excess ice thickness that exceeds the limit is computed as follows:

if (do_ice_limit .and. (sum(e2m) > max_ice_limit*DI*LI)) then
heat_limit_ice = sum(e2m)-LI*DI*max_ice_limit


do_ice_limit is a namelist logical parameter controlling whether to apply a limit on the sea ice
max_ice_limit is a namelist parameter, the specified limit [m], default of 4 meters.

The constraining heat (due to ice restoring or to an ice thickness limit) is then applied to the sea ice. For the simpler “slab ice” case, this is done as follows (more complicated treatment for non-slab ice is given in the ice model code):

tot_heat = heat_res_ice+heat_limit_ice
if (slab_ice) Ice%h_ice (i,j,2) = Ice%h_ice (i,j,2) - tot_heat/(DI*LI)

The restoring or ice limit heat fluxes are then archived as output from the model run, as a climatology of those fluxes normally will need to be computed for use in subsequent control or perturbation runs. The following equations convert the energy from (J/m^2) to flux (W/m^2).

Ice%qflx_lim_ice(i,j) = heat_limit_ice / dt_slow
Ice%qflx_res_ice(i,j) = heat_res_ice / dt_slow

5. Reconciliation of observed SST and sea ice climatologies:

The default “target” observations of SST and sea ice (concentration and thickness) that are used during the restoring run come from netcdf input data available as part of FMS. These were derived from the AMIP2 forcing data sets.

In practice, we wished to maintain a consistency between SST and sea ice observations for all times and locations in this target set of observations. To maintain this consistency even in the presence of spatial and temporal interpolation requires a consistency check be done following the space-time interpolations.

We let sea ice concentration be the primary variable in resolving conflicts between ice concentration, ice thickness, and sea surface temperature. Also, since the “slab ice”-type sea ice model allows only for zero or 100% ice coverage, a threshold must be used in deciding whether a grid point will be regarded as ice-covered or ice-free. We use 20% as the threshold concentration. At that concentration or above, the ice thickness is adjusted to be at least 1 meter thick. For the “slab ice” model, the ice concentration will be 100%, but for the dynamical sea ice model, the actual concentration, if at least 20% is retained for the model. Ice concentrations below 20% are set to zero with zero ice thickness. If a gridpoint is regarded as having sea ice (> 20% concentration) the SST is set to the freezing point of sea water. If a gridpoint is regarded as ice free, the SST is adjusted, if necessary, to be at least slightly above the freezing point of sea water.

The following code accomplishes these consistency checks. Note that although the concentration is not explicitly set to 100% for the “slab ice”-type ice model, elsewhere in the code if the “slab ice”-type ice model is used the ice concentration is assumed to be 100%.

where (icec >= 0.2)
    iceh = max(iceh, 1.0)
    sst = t_sw_freeze
else where
    icec = 0.0
    iceh = 0.0
end where
where (icec==0.0 .and. sst<=t_sw_freeze) sst = t_sw_freeze+1e-10

It should be noted that the namelist logical parameter mcm_ice for the routine ice_spec.f90 is set to “.false.” for the standard MCM slab ocean runs. That namelist parameter is used in coupled model work with the MCM specifically for fixed SST runs that mimic the old climate group GCM code (supersource).