The LES model, developed by Deardorff (1980) and Moeng (1984) for the atmospheric boundary layer, dynamically solves for unresolved (subgrid) turbulent kinetic energy (TKE). Several adaptations of this model are currently used to model oceanic boundary layers, including high-latitude deep convection in response to storm forcing (Harcourt et al 2001), and forced convection in more shallow free-slip ocean boundary layers (e.g. Skyllingstad and Denbo, 1995; Wang et al, 1996).
Following Moeng (1984), the advection terms of the momentum equations
are computed prognostically from vorticity components
Terms involving both subgrid and resolved TKE are absorbed into a
modified nonhydrostatic pressure 
. b is the buoyancy force, and the mean
acceleration is removed from the right-hand side of the vertical momentum
equation to maintain a prescribed mean vertical velocity.
The LES uses a finite series expansion of the equation of state about a
potential reference state 
with density 
chosen to
minimize truncation errors. After removing the horizontally averaged vertical
acceleration 
, buoyancy fluctuations 
account for
mean scalar dependence in the expansion coefficients, thermobaricity, cabbeling
and a ‘halobaric’ effect through 
, 
, 
and 
, as
modifications to the standard approximation of constant thermal and haline
expansion coefficients 
and 
. The salinity-related components of cabbeling and the
dependence of the haline expansion coefficient on mean salinity are small and
generally omitted.
The evolution of potential temperature is computed prognostically from
with corresponding equations for salinity and Lagrangian tracers. The velocity
components are subject to the condition of incompressibility, which is applied
to the momentum equations to solve diagnostically for 
. Unresolved
stresses are modeled in terms of the symmetric resolved strain rate and local
nonlinear eddy viscosity 
; scalar fluxes in terms of resolved gradients and eddy
diffusivity 
:
Parameterizations eddy viscosity diffusivity, 
and 
, depend upon subgrid TKE 
, its associated
dissipation length scale 
, turbulent eddy Prandtl number 
and constant
. The
length scale 
of unresolved motion is equal to the resolution scale 
, unless reduced
below that by stability or by proximity to a wall boundary. Parameterizations
for the dissipation length scale 
for stably stratified turbulence are discussed
below in the description of simulations. Resolution scale 
is governed
either explicitly by a spatial filter, implicitly by the numerical scheme, or
by a combination of the two.
The subgrid TKE equation
includes unresolved buoyant production 
as well as
unresolved shear production and transport, and parameterizes dissipation as
. When a
portion of the turbulent inertial subrange is resolved, constants 
and 
of the subgrid parameterization are determined by integrals
over the 
Kolmogorov energy spectrum for unresolved TKE smaller than
The thermodynamic ice-ocean interaction parameterizations used by McPhee
(2000) in 1-D Local Turbulence Closure (LTC) simulations of the Weddell Sea
Mixed Layer have been implemented in the LES to simulate the freezing and
melting of ice. A prescribed mean surface stress 
derived from
observations is applied at the upper model domain, and the velocity boundary
condition is nonslip in fluctuations with respect to the mean horizontal
velocity at the surface. Heat flux through the ice is also derived from ANZFLUX
measurements, and the model ocean heat flux at the upper boundary is
where the freezing point 
is the boundary temperature 
is the LES
temperature in the first grid level adjacent to the boundary, and 
is a
nondimensional Stanton number, adjusted for consistency with vertical grid size
.
Freezing and melting of ice is then determined by the imbalance between heat
fluxes in the ice and in the ocean.
Several LES runs have been carried out to test the stability of the upper ocean under conditions observed during ANZFLUX at certain times when ocean observations were not made. In general, they agree with the LTC results of McPhee (2000), showing the upper ocean becoming unstable to a process of thermobaric parcel detrainment, but at an earlier point in time. These simulations do not yet account for cabbeling effects, as this term was omitted in the preliminary simulations for comparison with LTC results.
Figures 1-3 below summarize some results from the most recent
simulation, for which a short animation is available online at
http://www.oc.nps.edu/thermobaricity under Numerical Experiments,
Simulation 3. This simulation was initialized from averaged typical profiles
measured in 1994 on year days 216 and 217, and forced with time-dependent
surface stress deduced from observations. High wind in the first two days
results in shear production of mixed layer TKE with classical forced mixed
layer scaling, 
(Figure 1, top). Thermobarically unstable parcels are
subsequently detrained from the mixed layer in plumes with with small vertical
TKE, large negative skewness (Figure 1, second from top) and positive values of
below
. Large
negative values of buoyant TKE production (Figure 1, third plot) resulting from
the conversion of TKE to entrainment work, are followed by significant release
of potential energy below the mixed layer as the detrained thermobaric plumes
become more negatively buoyant with depth. The fourth plot in Figure 1 tracks
the net buoyant production of TKE above and below 
separately, and
demonstrates that without thermobaricity the detrained plumes would be a net
sink, rather than a source of TKE. This sensitivity (lower right) to the second
order terms in the expanded equation of state arises from the high level of
compensation between the heat and salinity fluxes (lower left).
Including the cabbeling term in the equation of state will probably
delay thermobaric instability, owing to the negative sign of 
for the plumes.
Mean TKE below the mixed layer is relatively small, but the skewness of the
thermobaric plumes is large and negative, and individual plumes (Figure 2) have
large downward velocities concentrated in structures about 100m in horizontal
scale.
We propose that this preliminary work on LES model development be extended through quantitative verification comparisons with existing data sets from polar boundary layers below ice, and subsequently be used to make predictions about upper ocean instabilities for use in designing future experiments in the Weddell Sea.
Figure 1: The top figure shows vertical TKE 
integrated over
depth, separately within and below the mixed layer, along with a representative
one-third of the dimensional surface stress forcing scale 
. The profile
timeseries of mean skewness 
is illustrated in the second plot, and buoyancy
flux 
-the work rate for buoyant TKE production- is shown in the
third plot. The fourth plot shows the net work rate due to buoyancy production
of TKE separately within and below the mixed layer, both with and without
thermobaricity. The lower right profile also shows the effect of thermobaricity
on buoyant TKE production, averaged over the entire simulation, and the lower
left plot decomposes 
into heat and salinity flux components.
Figure 2: Vertical cross sections of potential temperature (left plot) and vertical velocity (right plot), on about day 220 of the simulation, illustrating a detraining thermobaric plume.
Preliminary LES model development has been guided by earlier development of the 1-D LTC model. While this has allowed some preliminary simulations, comparisons between the modified LES and observations of turbulence below the ice are necessary to verify the thermodynamic parameterizations of ice-ocean interaction. Data for verification will come from past observations of both ANZFLUX and SHEBA experiments. Comparisons will include vertical fluxes and dissipation rates as well as the evolution hydrographic profiles. Figure 3 demonstrates the equivalence of LES subgrid TKE dissipation and dissipation observations from a Turbulence Instrument Cluster (TIC) during the ANZFLUX experiment (n.b. two different time frames) by Tim Stanton and Miles McPhee.
Figure 3: Time series of measured dissipation profiles (left plot) from the two final days of observation East of Maud Rise, along with LES subgrid dissipation (right plot) from the simulation of the subsequent two days, initialized from observed hydrographic profiles. Model resolution is 6.25m.
The two quantities are equivalent above the instrument noise level because the LES resolution scale is the same order as the length scale characterizing the TIC response, and because most of the dissipation is at much smaller length scales. We anticipate three verification simulations, two based on ANZFLUX observations and one based on SHEBA results, each about one week in duration. Ice-ocean parameterizations will be tuned to give the best agreement in these comparisons.
Modeling for predictions of upper ocean instabilities will focus on what may have happened in the area East of Maud Rise after the end of ANZFLUX ocean observations. This region is of interest because hydrographic analysis indicates a weak barrier to thermobaric instability if the upper ocean layer, and because satellite observations indicate significant areas of open water (Drinkwater, 1997) in 1994, and large polynyas in other years. The goal of this modeling is to use the dynamics of the LES predict more accurately the point at which instabilities occur and the form they will take. Furthermore, timeseries from ensembles of virtual measurements embedded in the LES will predict the shape of observed thermobaric plumes and other features of the upper ocean, and assess the random and systematic uncertainties associated with different observing systems. Embedded measurements will include Lagrangian, isobaric, and tethered Lagrangian drifter ensembles as well as Eulerian point and volume-averaged timeseries representing ice-fixed measurements.
In short, these embedded measurements will be used to assess the effectiveness of an array of experimental strategies, thus aiding in the design of future field expeditions. These virtual measurement ensembles will simulate observations in about six LES cases of one to two weeks in duration each, initialized from different 1994 ANZFLUX CTD profiles from East of Maud Rise, and forced with observed time-dependent surface forcing also derived from observations.
The LES code is currently optimized for parallelvector machines, where up to twothirds of the load is scalable. Each LES case, with 105 - 106 grid points model up to several weeks of simulation time (about 105 timesteps), and will require approximately 800 cpuhours on a Cray J90 or SV1, or 200 GAU in NCAR's general accounting units. In core memory requirements are moderate at 30MWords per simulation. These computer resources will be obtained through applications for the use of both NCAR and DOD high performance computing resources that are available to carry out NSF funded research. Large Eddy Simulations will be carried out by Harcourt and Lherminier, LTC simulations by McPhee. Analysis of model output for verification and prediction will be done by Harcourt, Lherminier, Stanton and McPhee.
Deardorff, J. W. (1980). Stratocumuluscapped mixed layers derived from a threedimensional model. Bound. Layer Meteor. 18, 495-527.
Drinkwater, M. R. (1997). Satellite microwave radar observations of climate-related sea-ice anomalies. Bull. Am. Met. Soc., Proc. Workshop on Polar Processes in Global Climate. 13-15 Nov. 1996, 115-118.
Harcourt, R. R., E. L. Steffen, R. W. Garwood, and E. A. D'Asaro (2001). Fully Lagrangian floats in Labrador Sea deep convection: Comparison of numerical and experimental results. J. Phys. Oceanogr. in press.
Moeng, C.H. (1984). A largeeddy simulation model for the study of planetary boundarylayer turbulence. J. Atmos. Sci. 41, 2052--2062.
Skyllingstad, E. and D. Denbo (1995). An ocean largeeddy simulation of Langmuir circulations and convection in the surface mixed layer. J. Geophys. Res. 100, 8501-8522.
Wang, D., W. G. Large, and J. C. McWilliams (1996). Largeeddy simulation of the equatorial ocean boundary layer: Diurnal cycling, eddy viscosity, and horizontal rotation. J. Geophys. Res. 101(C2), 3649-3662.