GFDL - Geophysical Fluid Dynamics Laboratory

FMS Slab Ocean Model Technical Documentation

T. Knutson

Sept. 4, 2003

Revised May 20, 2009

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:

2. Background on Slab Model Experiments in Climate Studies

3. Procedure for Performing a 2xCO2 Sensitivity Experiment

     a. Restoring Run


b. Q-flux adjusted Control run


c. 2xCO2 perturbation run

4. Description of the Model Equations

5. 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

3. Procedure for Performing a 2xCO2 Sensitivity Experiment

In this section, the procedure for performing a prototypical climate sensitivity
experiment (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

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).