An operational ocean forecast system incorporating NEMO
and SST data assimilation for the tidally driven European
North-West shelf
Author List
O’Dea EJ
1
, Arnold AK
1
, Edwards KP
1
,Furner R
1
, Holt JT
2
, Hyder P
1
, Liu H
2
Martin MJ
1
,
Siddorn JR
1
, Storkey D
1
, While J
1
Affiliation
1
Met Office, Exeter, UK
2
National Oceanography Centre, Liverpool, UK
Abstract
A new operational ocean forecast system, the Atlantic Margin Model implementation of the
Forecast Ocean Assimilation Model (FOAM-AMM), has been developed for the European
North West Shelf (NWS). An overview of the system is presented including shelf specific
developments of the physical model, the Nucleus for European Modeling of the Ocean
(NEMO), and the Sea Surface Temperature (SST) data assimilation scheme. Initial validation
is presented of the tides and model SST. The SST skill of the system is significantly improved
by the data assimilation scheme. Finally, an analysis of the seasonal tidal mixing fronts shows
that these in general agree well with observation, but data assimilation does not significantly
alter their positions.
Lead Author’s Biography
In 2004, after completing a PhD. in computational fluid dynamics, Enda O’Dea joined the Met
Office to work in ocean modelling. He now develops ocean forecast models in the Ocean
Forecasting Research and Development (OFRD) group. Recently, the group has overseen
the transition from a POLCOMS based forecast system to a NEMO based forecast system for
the shelf seas around the U.K. Enda’s principal research area is in shelf seas forecasting and
interests include the dynamics of tides, seasonal stratification, shelf slope currents and
regions of fresh water influence.
Introduction
The North-West European shelf seas have been the subject of numerous
hydrodynamic models of increasing complexity and sophistication.
1,2
3D baroclinic
hydrodynamic models have evolved from research tools into operational forecast systems at
operational centres.
3,4
Recently, in addition to operational forecasts for hydrodynamic
variables, ecosystem models have been implemented operationally.
5
Such ecosystem models
provide model estimates of biogeochemical variables such as chlorophyll in complement to
remote Earth observations.
6
In global ocean and basin scale modeling, data assimilation has proved an invaluable
component for operational forecasting.
7
For the shelf seas however
8
, the necessary inclusion
of shorter temporal and spatial scale processes, in particular in relation to the interaction of
the tides and the shelf, has discouraged the widespread use of data assimilation in
operational systems. In this paper we outline the development of an operational modeling
system for both physical and biogeochemical parameters in the North-West European
continental Shelf (NWS) that includes assimilation of Sea Surface Temperature (SST) data.
In the Met Office’s Forecasting Ocean Assimilation Model (FOAM)
9
system, the core
dynamical model has recently
10
migrated to the Nucleus for European Modeling of the Ocean
(NEMO).
11
Using NEMO allows short-term operational ocean forecasting systems to employ
the same fundamental ocean model code as in the global and basin scale seasonal and
climate prediction systems at the Met Office. Adopting NEMO is also beneficial to the group
by becoming part of a large and active cross-institutional developer base. Like the open
ocean, the strategy for the assimilative shelf seas forecasting system described here is to
apply NEMO as the physics engine. Assimilation of SST adapts the existing open ocean
FOAM system in a manner suitable for application in the shelf seas. Whilst NEMO is also
coupled with the European Regional Seas Ecosystem Model (ERSEM)
12,13
for ecosystem
modeling, this paper details only the physics and assimilation. The coupled physics-
ecosystem will be the subject of a following paper.
13
Previously, operational modeling of the NWS at the Met Office utilised the Proudman
Oceanographic Laboratory (now National Oceanography Centre) Coastal-Ocean Model
System (POLCOMS)
14, 15
. NEMO was developed as an open ocean model and therefore
lacked many of the shelf-specific features found in models such as POLCOMS. It was
therefore necessary to incorporate significant modifications to NEMO to make it suitable to
replace POLCOMS as the operational model for shelf applications. It should be noted that the
POLCOMS system is a well established and validated system
16
and provides a reliable
reference system from which to compare forecast skill in any new modelling system.
The operational shelf seas forecasting system is run in the Met Office operational
suite on a daily cycle and forms part of the Europe-wide operational oceanography
contribution to Global Monitoring for Environment and Security (GMES). The NEMO
forecasting system documented in this work is now providing operational forecasts. Analyses
and five day forecast products for the NWS are provided as part of the MyOcean project, the
EC FP7 project that currently delivers the GMES Marine Core Service.
17
The remainder of this paper is divided into the following sections. First the model
domain and configuration are described. Second an overview of the model in general and
specific enhancements that have been developed to tackle the shelf seas dynamics in
particular are given. Third, the developments required for assimilation of SST on the shelf are
described. Fourth, a preliminary validation of the system is detailed. Finally, conclusions and
future developments are discussed.
SYSTEM DESCRIPTION
Physical context: The North-West European shelf
The Atlantic Margin Model (AMM) region shown, in Fig 1, covers the North-West
European shelf and part of the North-East Atlantic ocean. A key feature dividing the shelf from
the deep ocean is the shelf slope, running from Portugal to Norway. Associated with the shelf
slope is the important Joint Effect of Baroclinity and Bottom Relief (JEBAR)
18
process, which
drives a poleward shelf slope current common to many eastern margins. Examples of other
poleward eastern boundary currents in other regions include the coastal undercurrents of
Chile and California, the Alaskan slope current and the bottom layer shelf break current of
southwest Africa.
19
The shelf slope itself varies in width and steepness. It is particularly steep
along the Iberian slope to the west of Portugal and the Cantabrian slope to the north of Spain.
The combination of step bathymetry and terrain following (sigma) coordinates requires special
treatment for the modelling of Horizontal Pressure Gradients (HPG). The slope current is also
variable and along the Iberian and Cantabrian slopes it is seasonal. It manifests itself as the
Iberian Slope current during autumn and winter.
20,21
Further north the shelf widens from the Aquitaine slope to the Armorican and Celtic
shelf slopes with France to the east. The Celtic shelf slope is a major source of internal tides
22
and enhanced mixing.
23
Travelling northwards, the shelf slope encompasses the Celtic Seas
24
of the English Channel, Irish Sea, and the Celtic, Irish, Malin and Hebrides shelves.
Thereafter the shelf slope turns more eastward at the Faroes-Shetland ridge towards Norway
with the North Sea to the south and the Faroese channels to the north. Finally, crossing the
Norwegian trench to the south, the shelf slope travels parallel to the Norwegian coast towards
the northern boundary of the domain.
The Faroese channels and the Wyville-Thomson Ridge are important areas for the
return of the cold dense outflow from the Nordic seas.
25,26
Accurate modelling of the overflow
requires higher resolution models
than is currently computationally possible for operational
coupled physical–ecosystem models of the entire region. Careful attention is required at such
overflows particularly in relation to the limitations of terrain following coordinates and steep
bathymetry.
Although the main aim of the shelf model is in simulating the on-shelf properties, the
off-shelf dynamics and the shelf slope current are also important as they impact cross-slope
transport. Approximately 12.5% of the global tidal energy is transmitted into the Celtic Seas
from the North Atlantic,
with large tidal responses in the English Channel, Bristol channel and
Irish Sea.
24
The large tidal response results in large dissipation of tidal energy and an input of
turbulent kinetic energy into the water column. The seasonal variation in both wind, which
further adds to the production of turbulence, and heating, which adds buoyancy, leads to
seasonally stratified and mixed regions. The balance between mixing created by wind and
tide, and stratification by thermal heating, leads to tidal mixing fronts
at the boundary between
well mixed and stratified water columns.
27
Thus accurate representation of tidal dynamics,
turbulence production and dissipation, and the air-sea flux of momentum and heat are critical
for modeling the regional dynamics. The effects of winds are not limited to turbulence
production but also drive currents, which along with the buoyancy field provide the residual
circulation.
28
Furthermore, wind forcing in combination with atmospheric pressure can
produce large and potentially dangerous storm surges in the North Sea.
29
As such, a regional
model must include the atmospheric pressure gradient forcing, and the interaction of tides
and surges.
Sources of freshwater influence the baroclinic flow with inputs from rivers such as the
Rhine leading to dynamically complex Regions Of Freshwater Influence (ROFI)
30
and coastal
currents. Sources of low saline water are not restricted to local riverine sources alone. Low
saline water of Baltic origin, which may be considered as a large estuarine source, exchanges
with relatively high saline North Sea water flowing into the Baltic in a dynamically complex
transition area. The connection with the Baltic consists of the shallow sills and narrow straits
of the Kattegat, the Sound, the Great Belt and the Little Belt. To resolve the flow through
these channels requires relatively high resolution models
31
; where the resolution of the shelf
model is not adequate to resolve them, fluxes between the Baltic and Kattegat must be
specified as a special form of boundary condition. The dynamics are further complicated by
the intrusion of the relatively deep Norwegian trench as far as the Skagerrak. This guides
North Atlantic water along its slope into the Skagerrak, where it upwells and re-circulates.
32
The outgoing current flowing along the Norwegian coast consists of both low saline Baltic and
coastal water, and mixed North Sea and North Atlantic water. Thus any model of the region
must be able to represent the complex combination of haline, bathymetric, heating, tidal and
surge effects that all interplay in this region.
Physical Model
The model is designed to provide simulations of the on-shelf hydrodynamics,
biogeochemistry and light environments of the NWS. The high socioeconomic interest in the
area has led to an intensive modeling effort, with a variety of high-resolution models exploring
specific dynamical regimes in detail. However, in the context of an operational forecast
system that is coupled to sediment and ecosystem models, a regional approach that
interconnects the variety of dynamical regimes is required.
15
The existing coupled
POLCOMS-ERSEM Medium Resolution Continental Shelf (MRCS)
5
system is nested into the
physics only 12km POLCOMS-AMM model and has a resolution of approximately 7 km. The
new FOAM-AMM system extends the coupling of ERSEM outwards from the MRCS domain
to cover the entire AMM region from 40˚S, 20˚W to 65˚N, 13˚E. Thus FOAM-AMM replaces
both the existing POLCOMS operational models, POLCOMS-AMM
33
and MRCS with a single
domain.
Fig 1 depicts the North-West Shelf Operational Oceanographic System (NOOS)
bathymetry covering the AMM region, which is a combination of GEBCO 1’ data and a variety
of local data sources from the NOOS partners. The shelf break in Fig 1 has been highlighted
by the 200 m isobath. Also depicted is the perimeter of the existing MRCS
16
domain. In order
to ensure that the cross-slope exchanges of momentum and tracers are well represented a
hybrid s-
σ
34
terrain following coordinate system is employed in those models in order to
retain vertical resolution on the shelf, while allowing a reasonable representation of deep
water processes.
The resolution of FOAM-AMM is 1/15° latitude by 1/9° longitude. The horizontal
resolution of ~7 km lies between typical shelf wide resolution (~12 km) and high-resolution
limited-area models (1.8 km) sufficient to resolve the dominant fine scale physics on the
shelf.
28
It is not sufficient to resolve the internal Rossby radius on the shelf, which is of the
order 4 km, but well resolves the external Radius (~200 km). Ideally the model would be of
sufficient resolution to resolve both the internal and external radii, i.e. a resolution of the order
<2 km. At present the computational cost of such a system make this impractical for coupled
hydrodynamic-ecosystem operational forecasts.
The model bathymetry of POLCOMS is derived from the NOOS bathymetry. Some
smoothing was applied to steep bathymetry such as the shelf break in the derivation of the
existing 12 km AMM domain. This was to reduce HPG errors and improve the shelf slope
current.
33
The 12 km POLCOMS-AMM bathymetry has been interpolated onto the
replacement 7 km FOAM-AMM grid for inter-comparisons of the two systems. This ensures
that the new physics does not gain advantage simply by having more bathymetric information
in initial hindcast comparisons. In the initial validation stages NEMO-AMM was also run on the
12 km AMM grid for complete like for like comparisons at equivalent grid resolution.
The version of NEMO used in FOAM-AMM is v3.2. As it is necessary to model tides
and surges, a non-linear free surface is implemented using a variable volume
35
and time
splitting methodology, using ‘leap-frog’ time stepping. The corresponding baroclinic time step
is 150 seconds and the barotropic sub-cycle time step is 5 seconds. The momentum
advection is both energy and enstrophy conserving.
36
The lateral boundary condition on the
momentum is free-slip. Horizontal diffusion of momentum is specified using both Laplacian
and bilaplacian operators. Because FOAM-AMM utilises terrain following coordinates, it is
necessary for the specification of Laplacian diffusion to be applied on geopotential surfaces to
prevent spurious mixing in the vertical, and bilaplacian diffusion to be done on model levels to
retain stability. The coefficients of laplacian and bilaplacian diffusion are 30.0 m
2
s
-1
and
1.0x10
10
m
4
s
-1
respectively. The total variation diminishing (TVD) scheme is used for tracer
advection
37
. The tracer diffusion operator is only Laplacian and operates along geopotential
levels. The tracer diffusion coefficient is 50 m
2
s
-1
.
There are 32 terrain following coordinates in the vertical. The terrain following
coordinate system is modified in two important ways. Firstly, as in
34
and
3
, coordinates
transition from a stretched S-coordinate system in the deep to a uniform σ-coordinate system
on the shelf. Following
33
the critical depth
c
h
is defined at 150 m and the stretching
parameters are defined as θ=6 and B=0.8. Focused resolution in deep water at the surface is
important for air sea fluxes of heat, freshwater and momentum, and the bottom in relation to
the bottom boundary layer and bed friction. On the shelf uniform coordinates are preferred, as
in shallow regions very small vertical cells will tend to result in violations of the vertical
Courant–Friedrichs–Lewy (CFL)
38
condition.
An additional modification to the coordinate system is based upon a z*-σ approach.
39
A
major constraint on terrain following coordinates occurs when adjacent ocean depths differ
significantly leading to errors in the calculation of the HPG term.
40
To reduce the error the
initial S-σ-system is created using a smoothed envelope bathymetry rather than the input
bathymetry itself. The motivation of the smoothing is to limit the steepness of the model levels
to a given threshold. The threshold in FOAM-NEMO is chosen as 0.3. Thus for any two
adjacent depths,
env
ji
h
,
!"
env
ji
h
1, +
in the envelope bathymetry, the relative difference is:
3.0
,1,
,1,
<
+
−
=
+
+
env
ji
env
ji
env
ji
env
ji
hh
hh
rn
(1)
The smoother only deepens, it does not shallow, the envelope bathymetry relative to
the source bathymetry. The S-σ-coordinate system is then created based on the envelope
bathymetry. However, h is then masked for any grid cell that is lower than the input
bathymetry. Hence, the coordinate slopes are never more than a desired threshold, at the
expense of some vertical levels near steep bathymetry. At such points the levels intersect the
bed and levels are lost. Fig 2 depicts a z*-σ -coordinate system with its underlying smoothed
envelope bathymetry. A hybrid z*-S-σ system provides one way of reducing HPG errors, that
has the distinct advantage that the shape of the topography is not overly distorted by
bathymetric smoothing. However, the underlying HPG scheme must also be suitably posed to
minimize spurious velocities and cross-pycnocline mixing. Such errors result from inclined
model surfaces relative to both geopotential surfaces and isopycnal surfaces. Furthermore,
the z*-S-σ system does not resolve the issue of having a non-uniform surface box, which has
implications for surface fluxes.
The standard HPG schemes in NEMO were found to give unacceptably large errors
with the non-linear free surface, generating large erroneous velocities over steep topography
such as the shelf break. Furthermore, these schemes were not able to deal with the hybrid
vertical coordinate. To address this, a new HPG scheme was developed employing a
pressure Jacobian method rather than the widely used density Jacobian method. This can be
illustrated with the following formula:
!P
!x
z
=
!P
!x
S
"
!P
!z
!z
!x
S
(2)
Here, P is the pressure, z is the non-transformed physical vertical coordinate and s is the
transformed vertical coordinate used in the model.
z
x∂∂ /
refers to the partial derivative in the
horizontal defined on geopotential surfaces and
S
x∂∂ /
is the horizontal partial derivative
defined on coordinate surfaces of the model coordinate system. A constrained cubic spline
(CCS) method has been employed here to reconstruct the vertical density profile. The CCS
reconstruction has the property of monotonicity. The vertical pressure profile can be
calculated analytically, so the density Jacobian method is not needed. By splitting the second
term of the two-term pressure gradient formula into left and right hand side parts, the pressure
gradient can be calculated on the velocity cells without any weight parameter. In this formula,
there is no hydrostatic consistency constraint. This pressure Jacobian HPG method can be
applied to any hybrid vertical coordinate. For details about this HPG method, please refer to
41
. The combination of the new HPG scheme and the vertical coordinate scheme give good
results in proximity to steep topography.
The non-linear free surface allows for the accurate representation of tides and surges.
At the open boundaries tidal energy enters the domain via a Flather
42
radiation boundary
condition. Fifteen tidal constituents, calculated from a tidal model of the North-East Atlantic
43
,
are specified for the depth mean velocities and sea surface elevation. As the AMM region
covers a significant area, the equilibrium tide is also specified. In addition to the tidal
boundaries, FOAM-AMM is one-way nested within the Met Office operational FOAM 1/12˚
deep ocean model for the North Atlantic.
10
Temperature and salinity are relaxed to the values
specified by FOAM 1/12˚ model over a ten point relaxation zone on the open boundaries
using the flow relaxation scheme.
44
Sea surface elevation and barotropic currents from the
FOAM 1/12˚ North Atlantic model are added to the tidal constituents via the Flather boundary
condition.
Vertical turbulent viscosities/diffusivities are calculated using the Generic Length Scale
(GLS) turbulence model.
45
This allows for a choice from a range of closure schemes. In
FOAM-AMM, the second-moment algebraic closure model of Canuto
46
is solved with the two
dynamical equations
47
for the turbulence kinetic energy (TKE), k, and TKE dissipation, ε.
48
