Model Description: The Numerical Solution

The Grid

The region of interest is divided into a latitude-longitude grid with $(i,j)=(1,1)$ being at the southwest corner. The variables are arranged on an Arakawa B-grid, so both components of velocity are defined (at $u$-points) half a grid spacing to the southwest of the points where scalar variables are defined ($b$-points). The open and land boundaries of the model lie along $b$-points to remove the need for a horizontal velocity boundary condition at the land boundaries. In later versions of the model, particularly those with wetting and drying, this arrangement of the land boundary is modified.

The water column at each grid point is divided into $N$ (typically the order of 20 or 30) $\sigma$-levels and Fig. 1 shows a schematic of the levels at a grid point. These are indexed $k=0..N-1$, with levels $k=0$ and $k=N-1$ lying respectively below the sea bed and above the sea surface to facilitate the use of flux boundary conditions, hence levels $k=1..N-2$ are interior to the model. To maintain resolution near the surface in deep water we allow the spacing of $\sigma$-levels to vary in the horizontal, according to the transformation

$$\sigma_k = S_k+\frac{h_{i,j}-h_{c}}{h_{i,j}}(C(S_k)-S_k) h_{i,j}>h_{c}$$

$$= S_k h_{i,j} \le h_{c},$$ (27)

where $S_k$ are $N$ evenly spaced levels between $\sigma=-(1+0.5(N-2)^{-1})$ and $\sigma=0.5(N-2)^{-1}$. The deviation from the usual $\sigma$-levels is given by

$$C(S_k)=(1-B)\frac{\sinh(\theta S_k)}{\sinh(\theta)}+ B \fra... ...h(\theta [S_k+0.5])- \tanh(0.5 \theta)}{2 \tanh(0.5 \theta)}.$$ (28)

This follows the s-coordinate transformation of Song and Haidvogel (1994), except here we retain the definition of $\sigma$, so the spacing of the levels in $\sigma$-space does not vary with time. In Holt and James (2000) the parameters $h_{c}=150$m, $\theta=5$ and $B=0.25$ were used to give an increased resolution only at the surface. For example at a depth of 1800m there are then 7 levels shallower than 200m, in contrast to only 2 when evenly spaced $\sigma$-coordinates are used. The equations of motion require no further transformation when these levels are used since the advection scheme employed (section 3.3) is a finite volume formulation and the pressure gradients are calculated by interpolating the pressure onto horizontal planes (section 49.1).

The Depth Mean Solutions

The equations for the surface elevation (Eqn. 14) and depth mean currents (Eqns. 4 and 5) are integrated forward in time (using centred space forward time differencing) with a barotropic time step which is a fraction of the (long) baroclinic time step $\Delta t$. The surface elevations are filtered using the method described by Killworth et al. (1991) to prevent the "checkerboard" pattern of grid-scale noise which is sometimes found when using the B-grid.

Advection

In order to maintain horizontal gradients and minimise numerical diffusion the advection of momentum and scalars uses the `Piecewise Parabolic Method' (Colella and Woodward 1984, James 1996). This scheme assumes the variables vary parabolically across the grid boxes, with the mean in the grid box, $\overline{a}$, taken to be the value at the grid-point, $a_{i,j,k}$ . The parabolas are then defined by

$$a( \eta ) = a_L + \eta (a_R-a_L+a_6(1-\eta)), \quad a_6=6 \overline{a}-3(a_L+a_R),$$ (29)

where the coordinate $\eta$ varies from 0 to 1 across the grid box. The left and right values, $a_L$ and $a_R$, are found initially by interpolating a quartic polynomial fit to the integral of $a$. These are adjusted so that no new extrema are introduced by the parabolas. The advective flux is then calculated by integrating the parabolas in an upwind sense. After horizontal advection, the PPM scheme is used in the vertical to adjust the $\sigma$-levels back to their original spacing; this determines the flux across the vertical faces of the grid boxes. For scalar variables the depth mean flux used in the advection step is the sum of the fluxes calculated during each of the barotropic time steps (including the equivalent flux for the Killworth filter). This ensures that the depth at the end of the advection step is identical to that at the end of the barotropic step, so scalars are conserved during advection.

Horizontal Pressure Gradients

A novel method of estimating the horizontal pressure gradient term is used since the model domain may include regions of extremely steep topography at the shelf edge. The traditional method of calculating the horizontal pressure gradient in $\sigma$-coordinate models is to estimate the pressure (or density) gradients along $\sigma$-levels and then correct these for the slope of the coordinates. The source of error in this is well documented (Haney 1991, Mellor et al. 1994): in the case of a flat thermocline overlying steep topography both these terms are of the same magnitude and should cancel to zero, but this generally leads to a significant truncation error which can drive erroneous currents. A number of methods of overcoming this problem have been suggested and successfully applied, for example by Stelling and van Kester (1994) and Slørdal (1997). The approach adopted here involves estimating the four pressure gradients at the edges of the horizontal plane (an example is shown in Fig. 1) with $(u,v)_{i,j,k}$ defined at the centre and the corners being at the $b$-columns $(i,j)$,$(i-1,j)$, $(i-1,j-1)$ and $(i,j-1)$; for simplicity these are labeled $q=1..4$, and the vertical index for the $u$-point is dropped. The depth of the plane is

$$z=\frac{1}{4} (\sigma \sum\limits_q H_q + \sum\limits_q\zeta_q )$$ (30)

and at $b$-column, $q$, this is at $\sigma_q =(z-\zeta_q)/H_q$. The nearest $\sigma$-levels, $k_q$, above the plane (shown as double circles in Fig. 1) lie at a fractional distance, $r=(\sigma_{k_q} - \sigma_q)/(\sigma_{k_q} - \sigma_{k_q-1})$, from the plane. If the buoyancy is taken to vary linearly between $\sigma$-levels, the pressure due to this interval can be estimated and the pressure at $\sigma_q$ can then be written as

$$\psi_q= -0.5 \delta \sigma [ b_{q,N-2} + \sum\limits_{k=k_q}^{N-2} (b_{q,k} +b_{q,k+1})$$

$$+(2r-r^2)b_{q,k_q}+r^2b_{q,k_q -1} ],$$ (31)

for $1 \le k_q \le N-2$, and

$$\psi_q= -0.5 \delta \sigma [(2r-r^2)b_{q,k_q}+r^2b_{q,k_q -1} ],$$ (32)

for $K_q=N-1$; in this calculation level $N-1$ is at the sea surface, where the buoyancy is defined by $b_{q,N-1}=b_{0 q,N-2}$.

The pressure gradients along the edges of the plane are then

$$\Delta \psi_1 = \psi_1-\psi_2, \quad, \Delta \psi_2 = \psi_2... ...a \psi_3 = \psi_4-\psi_3, \quad \Delta \psi_4 = \psi_1-\psi_4$$ (33)

and the corresponding changes in velocity are

$$\Delta u_{i,j,k} =-0.5 \frac{\Delta t}{Rcos \phi \Delta \chi} (\Delta \psi_1 +\Delta \psi_3)$$

$$\Delta v_{i,j,k} =-0.5 \frac{\Delta t}{R \Delta \phi} (\Delta \psi_2 +\Delta \psi_4)$$ (34)

This technique allows a straightforward treatment of cases where the sea bed lies above the plane at one or more of the surrounding $b$-points; if either of the corners on an edge are below the sea bed then the corresponding pressure gradient takes the last defined value above it. This assumes the isopycnals are horizontal close to the sea bed and does not reflect the correct boundary condition for temperature and salinity (zero normal gradients). However, the diffusive layers in which the thermoclines and haloclines curve to meet this condition are not well resolved in this model.

This method has been tested against the pressure gradient due to an analytically defined thermocline and gives significantly more accurate results than the conventional method when the thermocline is flat or sloping in the opposite sense to the $\sigma$-levels. Both methods give equally good results when the thermocline slopes in the same sense as the $\sigma$-levels.

Vertical Mixing

The vertical diffusion of scalars and momentum uses a fully implicit numerical scheme and in the case of momentum, the surface and bottom stresses are introduced as fluxes at this stage. The diffusivities used in these calculations are estimated from Eqn. 25, with $l$ limited by stratification (Eqn 22). The equation for $q^2$ (Eqn. 18) is integrated forward in time with the shear, buoyancy and dissipation terms being split into implicit and explicit parts by a Taylor expansion about the current time level following Annan (1999). The diffusion term is solved implicitly with values of $K_q$ estimated between the values of $q^2$ by arithmetic averaging.

Convective Adjustment

As suggested by Deleersnijder and Luyten (1994), we find that the turbulence closure described in section 2.5 provides an inadequate treatment of statically unstable layers. This is not unexpected since convection is not a sub-grid scale process, so is not well represented by diffusion at a grid-point. Moreover tiny values of $N^2<0$ are found to give very large diffusivities and under surface cooling conditions these propagate down the water column but take a number of time steps to mix the surface layer.

Instead of the turbulence closure scheme we use a simple convective adjustment procedure to mix unstable density gradients more effectively. At each horizontal grid point the model searches down in the vertical for a point where $N_{k-\frac{1}{2}} ^2<0$, at $\sigma$-level $k=k'$. The temperature and salinity at $k'$ is then repeatedly averaged with $\Delta k$ points below until the resulting buoyancy $b_0( \overline{T}, \overline{S})_{k'-\Delta k} > b_0(T,S)_{k'-\Delta k-1}$ (here an overbar indicates an average from $k'$ to $k'- \Delta k$). The temperature and salinity in this interval is then replaced by $\overline{T}$ and $\overline{S}$. If necessary this procedure is repeated for points above $k'$. The velocity field is also adjusted by convection; prior to the adjustment at b-points, the $T$ and $S$ values are interpolated onto the u-points and $(u,v)$ are averaged in the vertical according to the above scheme.

Because of its arbitrary nature, the use of this convective adjustment scheme is less than satisfactory. This is particularly the case for velocities; while its effects on their profiles are not great it may result in excessive mixing under unstratified conditions. However, we do find the use of this adjustment somewhat improves the comparison with SST observations, but a more satisfactory method of treating statically unstable conditions needs to be sought.

The Initial and Boundary Conditions

The model may be initialised from rest with a three-dimensional climatology or other data field of temperature and salinity. The elevations use a flux/radiation boundary condition, with the external elevation and depth mean velocity given: in several applications this has been determined from 15 tidal constituents taken from a tidal model of the northeast Atlantic. A surge or density component of the boundary currents and elevations may be included if sufficient information to specify these is available. The effect of large scale atmospheric pressure gradients (the `inverse barometer effect') is not included in the boundary conditions. The implications of omitting these boundary components are discussed in Holt et al. (2001).

Temperature and salinity may be relaxed to climatological values in a region next to the model boundaries. These values, $T0$, are taken from the boundary data. We find these temperature and salinity fields often contain some small scale variations, for example on the continental slope, which can be removed by smoothing on horizontal planes since these variations cannot adjust and may drive erroneous currents. The effects of this smoothing are examined in Holt et al. (2001). For a relaxation zone 4 points wide, the temperature at the grid points $p=0\ldots 3$ from the western boundary are given by

$$T_{i+p,j}=r T0_{i+p,j}+(1-r) T_{i+p,j},\quad r=(4-p)/4,$$ (35)

with similar equations for other boundaries and for salinity. The corresponding relaxation time scale is $t \sim \Delta t/r$; about 26 minutes for $p=3$.

The sea surface is forced by data from a weather prediction model and the temperature updated by the heat fluxes calculated from bulk formulae using this data following Gill (1982) and downwelling solar radiation (with a decay scale taken to be 0.154m$^{-1}$). Details of these are given by Holt and James (1999a) and will not be repeated here. At present there is no surface salinity flux in this model; however fresh water inputs are determined from daily discharge data from rivers plus, for applications including the North Sea, the exchange with the Baltic and the Kattegat.