POLCOMS technical description

Introduction

This section of the POLCOMS documentation is a general description of the model, including the model equations, the model forcing and the numerical methods employed. Much of this material has been published in Holt and James (2001). It ends with a bibliography including published work using POLCOMS.

The physical model used in POLCOMS is based on the Proudman Oceanographic Laboratory three-dimensional baroclinic model POL3DB and has been developed to incorporate features suitable for the modelling of baroclinic processes on the shelf, at the shelf-slope and in ocean regions to allow long term coupled ocean-shelf simulations. Versions of this model have been running operationally at the UK Met Office since 2000, though its development can be traced back to the mid 1980s, as shown in the bibliography.

The simulation of seasonal currents, water quality parameters and plankton variability in shelf seas requires a physical model which goes beyond the well established tide and surge models to include the effects of horizontal and vertical density variations. This is necessary since the former introduces seasonal transports not seen in constant density models, while the latter controls the vertical fluxes crucial to biological production. Moreover thermal fronts affect both the horizontal and vertical transport of tracers. To this end a three-dimensional model with temperature and salinity treated as prognostic variables has been developed: the Proudman Oceanographic Laboratory Three-dimensional Baroclinic B-grid model (POL3DB). The origins of this model lie with James (1986) and it has subsequently been developed (James 1996) to include a sophisticated advection scheme, the `Piecewise Parabolic Method' (PPM) described below. This has excellent feature-preserving properties making it ideal for the simulation of on-shelf baroclinic features such as river plumes (James 1997) and fronts (Proctor and James 1996), and the transport of tracers from localised sources (Holt and James 1999b). Moreover, the model is formulated on an Arakawa B-grid (see section 3.1), in contrast to the C-grid used in many shelf sea models, for example the Princeton Ocean Model (Blumberg and Mellor 1987). The B-grid is more commonly used in deep ocean models, for example Killworth et al. (1991) and is well suited to the modelling of horizontal density variations since the Coriolis term can be calculated without averaging. This helps prevent the dispersion of the velocity features associated with fronts, in contrast with the C-grid which requires averaging over a number of points to calculate this term. A disadvantage is that continuity and scalar advection require averaging of velocities not required on the C-grid. It is this choice of grid and the non-diffusive advection which particularly distinguish this model from others.

Since the model is designed to cover sea areas which include both shelf and deep-sea regions it includes a number of techniques pertinent to the treatment of deep water: although the model equations remain in $\sigma$-coordinates in the vertical, the spacing of these coordinate surfaces on the finite difference grid is allowed to vary in the horizontal according to the s-coordinate transform of Song and Haidvogel (1994) (see section 3.1); the pressure gradient calculations are made by interpolation onto horizontal planes through the points where velocities are defined (section 49.1); and a term for the variation of compressibility with temperature and salinity (section 2.3) is included in the model equation of state for sea water. Horizontal diffusion is included as an option, but it is not needed for model stability.

Model Description: The Model Equations.

The model is formulated in spherical polar sigma coordinates: $\chi$ (eastwards), $\phi$ (northwards) and $\sigma$ (vertical), with $\sigma=(z-\zeta)/(h+\zeta)$, where $z$ is the Cartesian vertical coordinate, h is the water depth relative to the reference sea level ($z=0$) and $\zeta$ is the elevation above this; the total water depth is H=h+ $\zeta$. While $\sigma$ is the vertical coordinate it is discretised onto levels which vary in the horizontal in $\sigma$-space; this is described in section 3.1.

The Equation of State

The density is defined by an approximation to the full UNESCO equation of state: $\rho(T,S,p)=\rho(T,S,0)+\rho'(T,S,p)$, where $T$ is the potential temperature ($^\circ$C), $S$ the salinity (p.s.u.) and $p$ the pressure relative to the sea surface. $\rho(T,S,0)$ is taken from the UNESCO equation of state and

$$\rho'(T,S,p) = 10^4 \frac{p}{c^2} (1-0.20 \frac{p}{c^2}),$$ (1)

with

$$c = 1449.2+1.34(S-35)+4.55T-0.045T^2+0.00821p+15.0 \times 10^{-9} p^{2},$$ (2)

following Mellor (1991).

For accuracy in the numerical calculation we define the buoyancy $b=b_0+b'$, where the 'potential' buoyancy is $b_0=g(\rho_0-\rho(T,S,0))/ \rho_0$ ($\rho_0=1027$kgm$^{-3}$ is the reference density) and the variation of compressibility with temperature and salinity is accounted for by $b'=g(\overline{\rho'(Z)}- \rho')/ \rho_0$, with $Z=z-\zeta=\sigma H$. Our initial condition gives $\overline{\rho'(Z)}=-0.004564Z$. The total (hydrostatic) pressure is then given by

$$P= P_a+\rho_0(\psi + g \zeta - gz)+0.002282gZ^2,$$ (3)

where $P_a$ is the atmospheric pressure and $\psi= H \int_0 ^{\sigma} b d \sigma$.

The Equations of Motion.

We solve the incompressible, hydrostatic, Boussinesq equations of motion and to allow time splitting between barotropic and baroclinic components, these equations are divided into depth varying and depth independent parts; so the eastwards velocity is $u=\overline{u}(\chi,\phi,t)+u_r(\chi,\phi,\sigma,t)$ and the northward velocity is $v=\overline{v}+v_r$. The depth mean equations are

$$\frac{ \partial \overline{u}}{ \partial t}= f \overline{v} -... ...ial P_a}{ \partial \chi} \right] +H^{-1}[F_S -F_B]+NLB_{\chi},$$ (4)

and

$$\frac{ \partial \overline{v}}{ \partial t}= -f \overline{u} ... ...ial P_a}{ \partial \phi} \right] +H^{-1}[G_S -G_B]+NLB_{\phi},$$ (5)

while the equations for the depth varying components are:

$$\frac{ \partial u_r}{\partial t}= -L(u)+f v_r + \frac{uv \tan \phi}{R} - \Pi _{\chi} + D(u) - H^{-1}[F_S -F_B]-NLB_{\chi},$$ (6)

and

$$\frac{ \partial v_r}{ \partial t}= -L(v)-f u_r - \frac{u^2 \tan \phi}{R} - \Pi _{\phi} + D(v) - H^{-1}[G_S -G_B]-NLB_{\phi},$$ (7)

where $R$ is the radius of the earth and the buoyancy terms are given by

$$\Pi_{\chi}=(R \cos \phi)^{-1} \left. \frac{\partial \psi}{\p... ...-1} \left. \frac{\partial \psi}{\partial \phi} \right\vert _z.$$ (8)

These have not been transformed onto $\sigma$-coordinates (unlike in previous applications) because of the particular method of solution used here (section 49.1). The depth means of the nonlinear and buoyancy terms are

$$NLB_{\chi}=\int_{-1}^0 \left[ -L(u) + \frac{uv \tan \phi}{R} - \Pi_{\chi} \right] d \sigma,$$ (9)

$$NLB_{\phi}=\int_{-1}^0 \left[ -L(v) + \frac{u^2 \tan \phi}{R} - \Pi_{\phi} \right] d \sigma.$$ (10)

The advection terms are given by

$$L(a)= \frac{u}{R \cos \phi} \frac{\partial a}{\partial \chi... ...{\partial \phi} + \Omega \frac{\partial a}{\partial \sigma},$$ (11)

with

$$\Omega=-\frac{\sigma}{H} \frac{\partial \zeta}{\partial t}- ... ...\left( H \cos \phi \int_0^{\sigma} v d \sigma \right) \right].$$ (12)

As is common practice in models of this type the vertical gradients of the stresses, $(\tau_{z \chi},\tau_{z \phi})/\rho$, have been replaced by a diffusion term:

$$D(a)=H^{-2} \frac{\partial}{\partial \sigma} \left( K_z \frac{\partial a}{\partial \sigma} \right),$$ (13)

where $K_z$ is the eddy viscosity defined in section 2.5. The equation for the free surface is

$$\frac{\partial \zeta}{\partial t} = -(R \cos \phi)^{-1} \le... ...partial }{\partial \phi} (H \cos \phi \overline {v} ) \right].$$ (14)

We use slip vertical boundary conditions: the components of surface and bottom stress and the corresponding friction coefficients are given by

$$(F_s,G_s)=c_s \frac{\rho_A}{\rho_0} (u_w,v_w) \sqrt{ u_w^2+v_w^2}, \quad c_s =0.63 + 0.66\sqrt{ u_w^2+v_w^2}$$ (15)

following Smith and Banke (1975), where $(u_w,v_w)$ is the wind velocity at 10m; and

$$(F_B,G_B)=c_B (u_B,v_B) \sqrt{ u_B^2+v_B^2}, \quad c_B=\left... ...t( \frac{\delta}{z_0} \right) \right]^{-2}, \quad c_B>0.005.$$ (16)

following Blumberg and Mellor (1987). The near bed velocity, $(u_B,v_B)$ is defined at a depth $\delta$ above the sea bed, the roughness length is taken to be $z_0=0.003$m and $\kappa=0.41$ is von Karman's constant.

The transport equation for temperature $T$ (and salinity, $S$) is

$$\frac{\partial T}{\partial t} = -L(T)+D(T),$$ (17)

which uses the eddy diffusivity $K_H$ (defined in the next section).

Turbulence Closure

To estimate the vertical eddy viscosity and diffusivity we use the Mellor-Yamada-Galperin level 2.5 turbulence closure scheme (Mellor and Yamada 1974, Galperin et al. 1988), with an algebraic mixing length. This results in a one equation, second moment closure lying somewhere between level 2 (steady state) and level 2.5 (two equation). We have included in this scheme a number of modifications which, while not necessarily optimal, were found to improve the on-shelf comparison between model and observation (both CTD and AVHRR), in a particular application (Holt and James 1999a).

The evolution of $q^2$ (twice the turbulent kinetic energy density) is given by

$$\frac{\partial q^2}{\partial t} = 2 K_z M^2-2(K_H- \alpha K _z) N^2 - \frac{2 q^3}{B_1 l}+D(q^2)$$ (18)

where

$$M^2= H^{-2} \left[ \left( \frac{\partial u}{\partial \sigma... ...ight], \quad N^2= H^{-1} \frac{\partial b_0}{\partial \sigma},$$ (19)

and $B_1=16.6$. The term in $\alpha K _z$ (with $\alpha$=0.7) is a simple representation of vertical mixing by long wavelength internal waves following Mellor (1989).

There are many choices of algebraic mixing length, $l$, see for example Xing and Davies (1996), and we have not evaluated all these options, but instead use the Bakhmetev scale since Proctor and James (1996) and Holt and James (1999a) show this gives reasonable results in the North Sea. We modify the length scale in deep water so that $l=\kappa H l_0(\sigma)$; where $h\le h_c$:

$$l_0= (1+ \sigma) (- \sigma)^{0.5},$$ (20)

and where $h>h_c$

$$l_0= \frac{h_c}{h} \left(1+\sigma \frac{h}{h_c} \right) \left(-\sigma \frac{h}{h_c} \right)^{0.5} \sigma \ge \frac{-h_c}{3 h}$$

$$= \frac{2 \sqrt{3}}{9} \frac{h_c}{h} \frac{-h_c}{3h} > \sigma > \frac{2h_c}{3h}-1$$

$$= (1+\sigma) \left(1- \frac{h}{h_c} (1+ \sigma) \right) ^{0.5} \sigma \le \frac{2h_c}{3h}-1.$$ (21)

This maintains the surface and near bed mixing length profiles at a water depth of $h_c$, (=150m for this work) into deeper water, with intermediate depths in the water column taking the maximum value at $h=h_c$. Defining the profile in this way simulates, in an arbitrary fashion, the eddy scale becoming independent of water depth in deep water and prevents $l$ increasing without limit with $h$. The eddy size in stratified water is limited according to the Ozmidov length scale (for dissipation, $\epsilon$), $L=( \epsilon / N^3)^{1/2}$. Hence, we impose a limit on $l$ proportional to this, following Galperin et al. (1988):

$$l^2<0.28 \frac{q^2}{N^2}.$$ (22)

The near bed boundary condition arises from a steady state balance between shear and dissipation:

$$q^2 \vert _ {\sigma=-1} = B_1 ^{2/3} \sqrt(F_B^2+G_B^2),$$ (23)

whereas at the surface a balance between diffusion and dissipation is used to represent the flux of turbulence from surface wave breaking. Following equation 27 of Craig and Banner (1994), this gives

$$q^2 \vert _ {\sigma=0} = \alpha_{wave}^{2/3} (3B_1/S_q) ^{1/3} \sqrt(F_S^2+G_S^2),$$ (24)

with $\alpha_{wave}=100$. For a given stress this leads to an increase in the value of the surface turbulence by a factor of about 20 above the sea bed boundary condition. This crude representation of surface wave effects leads to an improvement in the summer thermocline depth, but is flawed in a number of respects. First, no spatial and temporal variation of the wave energy factor, $\alpha_{wave}$, is included and second the vertical resolution is unable to resolve the narrow region in which this surface turbulence is dissipated so the effects are allowed to penetrate deeper than the work of Craig and Banner (1994) would suggest.

This system is closed by

$$K_z= S_M l q \quad K_H= S_H l q \quad K_q=S_q l q,$$ (25)

using the stability relations from Galperin et al. (1988):

$$S_M=\frac{0.3933-3.086G}{(1-34.68G)(1-6.127G)} \quad S_H=\frac{0.4939}{1.0-34.68G} \quad S_q =0.2$$ (26)

where $G=-l^2 N^2/ q^2$ and Eqn. 22 gives $G>-0.28$. $K_q$ is the coefficient for the vertical diffusion of $q^2$. In unstably stratified conditions, a convective adjustment scheme is used to mix vertically before the diffusivities are calculated; this is described in the next section.