A nonhydrostatic numerical model for density stratified flow and its applications

706

2008,20(6):706-712

A NONHYDROSTATIC NUMERICAL MODEL FOR DENSITY STRATIFIED FLOW AND ITS APPLICATIONS* LI Qun Polar Research Institute of China (PRIC), Shanghai 200136, China, E-mail: [email protected] XU Zhao-ting Physical oceanography Institute, Ocean University of China, Qingdao 266100, China (Received September 28, 2007, Revised December 25, 2007) Abstract: A modular numerical model was developed for simulating density-stratified flow in domains with irregular bottom topography. The model was designed for examining interactions between stratified flow and topography, e.g., tidally driven flow over two-dimensional sills or internal solitary waves propagating over a shoaling bed. The model was based on the non-hydrostatic vorticity-stream function equations for a continuously stratified fluid in a rotating frame. A self-adaptive grid was adopted in the vertical coordinate, the Alternative Direction Implicit (ADI) scheme was used for the time marching equations while the Poisson equation for stream-function was solved based on the Successive Over Relaxation (SOR) iteration with the Chebyshev acceleration. The numerical techniques were described and three applications of the model were presented. Key words: nonhydrostatic assumption, stratified flow, self-adaptive grid, nonlinear internal wave

1. Introduction  Numerical models have become an indispensable tool for ocean modeling. Such models are typically based on the hydrostatic approximation in view the argument that their horizontal length scales are larger than the vertical length scales. There are a wide variety of physical processes in oceans, and many of these processes are adequately modeled with the hydrostatic approximation. However, hydrostatic models cannot be used to simulate small-scale (i.e., scale comparable to or smaller than the water depth) physical processes, including such important cases as shear instabilities, high-frequency internal waves, and wave breaking[1]. These small nonhydrostatic physics process play an important role in understanding the energy cascade in the ocean, for example, internal solitary waves contribute to the physics that influence mixing in a density stratified system and have been previously

* Project supported by the National Science and Technology Supporting Plan (Grant No. 2006BAB18B03). Biography: LI Qun (1981- ), Male, Ph. D.

shown to be non-hydrostatic. When non-hydrostatic pressure is neglected in a model, the governing equations may miss a piece of the physics (dispersion effect)[2-4]. Despite the fact that the non-hydrostatic pressure is necessary for correctly modeling the physics of some small scale process, the high computational cost of solving the nonhydrostatic pressure has limited its use in large-scale systems. So far, these nonhydrostatic geophysical problems have been studied mostly using pure two-dimensional models[5-7]. In this article, a nonhydrostatic numerical model is developed based on the vorticity-stream function equation for the study of density stratified flow over irregular topography. The layout of this article is as follows. Section 2 introduces the model equations, numerical scheme and the parameterization method. Three applications are described and model results are presented in Section 3. Some concluding remarks are given in the final section. 2. Model description The model is based on a two-dimensional Cartesian coordinate system with the x-axis lies at the

707

N 2  x, z , t , u z2  x, z, t

rigid-lid surface and the z-axis vertically upward. The rotation effect is included and the bottom topography depends only on the x coordinate, i.e. H = H  x

Ri( x, z, t ) =

with H y = 0 . The system of equations (under the

ab and kb are dissipation parameters describing

Boussinesq approximation) describing the dynamics of a continuously stratified fluid in an f-plane is expressed in terms of the stream function by[8]

Zt + J Z ,\  fvz = g

Ux +  ah\ xx xx + Ua

 ah\ xz xz +  av\ zz zz +  ah\ xz xz , vt + J Z ,\ + f \ z =  ah vx x +  av vz z ,

Ut + J  U ,\ +

^k

v

Ua g

`,

k h = P 0 + P  'x

ª¬ U  U 0  z º¼ z

z

(1)

is velocity, U , Ua , U0  z are wave induced perturbation of the density, average density, and undisturbed density profile, respectively, f is the Coriolis parameter, ah , av are respectively the horizontal and vertical coefficients of turbulent viscosity, while kh , kv the horizontal and vertical coefficients of density diffusion.

N  z

is the

Brunt-Vaisala frequency and J  , is the Jacobian operator, which is defined as J  a, b = ax bz  az bx . The coefficients av and kv are determined by the Richardson-number-dependent turbulent parameterization given by Pacanowski and Philander[9]:

a0 ª¬1+ D Ri  x, z , t º¼

p

a0 ª¬1+ D Ri  x, z , t º¼

(2) where

frequency ( N  x, z, t =  g U z U 0 ). This parameterization for the vertical turbulent kinematic viscosity and diffusivity increase the coefficients av and kv in regions with small values of Ri ( x, z , t ) , and in regions with vertical density inversions the Richardson number is taken to be zero. The horizontal eddy viscosity ah and turbulent diffusivity kh are parameterized by [10]

ah = E 0 + E  'x

where \ is stream function, Z is vorticity, u, v, w

kv =

parameters, N  x, z, t is the local Brunt-Vaisala

N 2  z \ x =  kh U x x +

Z = \ xx +\ zz , u = \ z , w = \ x

av =

background turbulence, a0 , D and p are adjustable

p

+ ab ,

+ kb

2

2

wu , wx wu wx

(3)

where E 0 and P0 are homogeneous background values, E and P are adjustable parameters, and 'x is the horizontal grid size. We are interest in the baroclinic motions of the ocean, and only the baroclinic response of the ocean will be considered. Thus the model assumes a rigid-lid approximation and assumes zero tangent wind stress on the free surface at z = 0 . Thus, on the free surface z = 0 we have

\ = G , Z = 0 , U z = 0 , vz = 0

(4)

where G is equal to zero or other forcing forms corresponding to particular configuration. For the viscous case at the bottom z = H  x the kinematics no-slip condition and the absence of heat and salt fluxes are assumed:

\ = 0 , z = H  x , v = 0 , Un = 0

(5)

where n is the normal to the bottom. The boundary value of the vorticity Z at the bottom is estimated with the use of the equation Z = ’ 2\ , where the stream function \ is taken from the previous temporal step. On the sidelong vertical boundaries the

708

baroclinic wave is assumed to vanish. For the purpose of flexibility, the Flow Relaxation Scheme (FRS) is employed and all the field variables are relaxed to the barotropic background state. Meanwhile, the modular techniques of FORTRAN 95 are explored in order to make the model flexible and scalable. Equation (1) will be numerically solved after a transformation from z-coordinate to a V -coordinate. Such a transformation will transform the irregular model region to a rectangular one.

DF wF wF = +Ui Dt wt wxi

(7)

where U i represents the velocity of the basic flow in the xi direction. The first term on the right hand side

z

x1 = x , z1 =

³ N  s ds

³

0  H  x

0

N  s ds

(6)

This transformation takes into account the fluid stratification as well as the variability of the bottom topography. Such a replacement is similar to the V -transformation frequently used in the numerical models of the ocean circulation. However, the basic difference is the transformation in Eq.(6) gives a fine resolution not only in the shallow part of the ocean but also at the pycnocline layer[8]. The splitting-up method (ADI scheme) is applied for the finite difference approximation of the equations. The vorticity transport equation is integrated in time by splitting the temporal step into two semi-steps. The spatial derivatives are approximated by means of second-order central finite differences. At each temporal semi-step an implicit system of equations with a tridiagonal matrix is obtained that is solved using standard techniques. The stream-function is then computed from the vorticity by solving the Poisson equation Z = \ xx + \ zz in the transformed V -coordinate. For its solution a standard relaxation method (SOR iteration with the Chebyshev acceleration) is applied. Finally the density perturbation is computed by using the same methodology as that for the vorticity. 3. Applications 3.1 Case 1: Generation of internal wave by the interaction of tide with topography Tidal flow over an obstacle induces vertical currents and their time variations excite internal waves. To investigate the wave excitation properly, it is necessary to consider the total time-variation of the forcing to which fluid parcels are subjected. By denoting the forcing as F , the total time variation of the forcing DF / Dt is given in Eulerian variables as [11]

Fig.1 Internal waves generated by a strong semi-diurnal tide over sea sill and its disintegration ( kU 0 V f  1 )

709

expresses the fact that the local time variations of the forcing are capable of exciting internal waves when the forcing is given by an oscillating tidal flow, and this term excites only internal waves of tidal frequency. On the other hand, the second term, which represents the advection effect, shows that the spatial variations of forcing in the frame fixed to the ground can also excite internal wave. This is because there is the temporal variation of forcing in the frame moving with the fluid at speed U i in the presence of the spatial variation of forcing in the fixed frame.

Fig.2 Internal waves generated by a strong semi-diurnal tide over sea sill and its disintegration ( kU 0 / V f !! 1 )

The first numerical test simulates the generation of internal waves by the interaction of barotropic tide with sill topography. The numerical model uses a horizontal grid size of 100 m, a vertical grid size of 5 m and a time step of 5 s. Figure 1 presents a time series of the internal tide generation process ( kU 0 / V f | 0.78 ). At the flood phase, a dome-like elevation of the pycnocline is formed over the ridge, this initial dome-like elevation is subsequently destroyed and separates into two progressive waves propagating in opposite directions. Due to the advection of the tidal current, an asymmetrical structure of the two progressive waves is formed. After further evolution, an internal solitary wave train begins to form, which is a result of the balance of the nonlinearity and the dispersion. In Fig.2, unsteady lee wave generation process ( kU 0 / V f | 1.8 ) is presented. Due to the strong advection effect of the tidal flow, the asymmetrical structure is clearer than in the first case. Previous numerical and theoretical work revealed that the generation process can be reasonably modeled with hydrostatic type models, however, for a precise description of the further evolution process, the nonhydrostatic effect should be considered. 3.2 Case 2: Internal solitary wave shoaling on a slope Remote sensing and in situ observations have revealed the frequent occurrence of large amplitude internal solitary wave on the continental shelf[12,13]. This kind of waves experience strong transformation during the runup process over the continental shelf [14,15] . One method for modeling the internal solitary wave evolution over variable topography is based on the weakly nonlinear and weakly dispersion KdV type equations, which has been widely used. Even though the KdV type model could successfully model the kinematics characteristics of the long nonlinear internal wave, it can not be used for modeling the strong interaction of internal solitary wave with bottom topography, especially with breaking happening. The second test simulates the propagation of large amplitude internal wave across a slope with water depths from 260 m to 100 m. The experiment set up is chosen to allow comparison with the observations from ASIAEX [16]. The background stratification is given by a theoretical result which gives a good approximation to the observed data during the experiment. The latter section is sufficiently long to enable the initial input to transform into a nearly steady solitary wave. The numerical model is used to simulate the propagation of the solitary wave across the continental slope from 260 m to 100 m over 20 km using a range grid size of 10 m, a vertical grid size of 1 m and a time step of 1 s. Figure 3 shows a time series of an internal solitary wave with amplitude of 50 m run-up the

710

continental slope. The simulations show the broadening of the initial wave with a decreasing forward propagating slope and the appearance of waves of elevation behind the main forwardly propagating wave. The numerical predictions are compared with observations from ASIAEX[16]. Figure 4 compares the predictions of the model with observations at a similar depth using ADCP during ASIAEX. This comparison shows reasonable agreement for both the wave profile and the distribution of the horizontal velocity (see Fig.4).

Fig.4

Fig.3 Internal solitary wave shoaling on the continental slope of the northern South China Sea

Comparison of the wave profile and horizontal velocity field with observations

The balance of the nonlinear steeping and nonhydrostatic dispersion makes it is possible for the existence of internal solitary wave. Therefore, it can be expected that the nonhydrostatic pressure gradient play an important role during the evolution of internal solitary wave. Since the baroclinic pressure gradient is the primary force term in maintaining the internal wave dynamics, the effect of Nonhydrostatic Pressure Gradient (NPG)[2] is evaluated by comparing the nonhydrostatic pressure gradient to the Baroclinic Pressure Gradient (BPG). The maximum NPG and BPG at each time step are plotted in Fig.5, together with the ratio of the maximum NPG and the maximum BPG. The baroclinic pressure gradient is greater than the nonhydrostatic pressure gradient, which illustrates that the baroclinic forcing is indeed the driving force for the soliton evolution. The nonhydrostatic pressure gradient is important in driving the soliton evolution as well, since it is approximately 4% to 30% of the baroclinic pressure gradient. Between 100 min to 200 min, the NPG, the BPG and the NPG-BPG ratio keep

711

approximately constant. There is a rapid increase in the NPG, the BPG and the NPG-BPG ratio between 200 min and 280 min, and then they gradually decrease with oscillations.

jump and propagates downstream along the bed (labeled “bottom-trapped vortex” in Fig.6). This vortex core has lighter waters in its core, but is trapped along the bottom due to pressure gradients that are stronger than the buoyancy forcing, and this result is similar to that in the previous theoretical work [17] .

Fig.6 Vorticity generation and boundary layer separation

Fig.5

Top: Time variation of maximum BPG and NPG. Bottom: Time variation of ratio (in percentage) of the BPG and NPG

The change in the magnitude of NPG, BPG, and the NPG-BPG ratio is mainly due to soliton topography interaction. Compared to the ambient region, stronger NPG and BPG exist within the soliton region. The peak occurs at the shelf break when the soliton goes through dramatic evolution due to wave-topography interaction (see Fig.3). 3.3 Case 3: Steady flow over Gaussian topography In the third test, we examine the ability of the presented model to simulate the steady stratified flow over vary topography. We simulate a field-scale, idealized, steady flow over a Gaussian topography. This type of problem has been studied in many papers, and is common in oceanographic (and even atmospheric) settings where currents force a stratified profile past topographic features. The twodimensional domain is 256 m deep by 5000 m long, with a Gaussian sill centered at 2250 m from the left side. The Gaussian topography is described 2 2 by, He  x / 2V in which V 2 = 105 and the sill height is 196 m. This makes the shallowest spot above the sill 60 m deep. The domain is forced by a constant inflowing current from the left face at 0.25 m/s (which is in the range of tidal currents). The stably stratified initial and inflowing density profile is given by U = 1001  e0.0673 z , where z = 0 is at the top of the domain. Figures 6 and 7 presents the typical numerical results of the vorticity field. A large vortex patch is shed from the shear layer (see Fig.6, the primary vorticity generation zone) at the internal hydraulic

Fig.7 Steady stratified flow past a Gaussian topography. Vorticity plots with a time interval of 750 s

4. Discussion and conclusion Based on the vorticity-stream function equations in the vertical two-dimensional coordinate systems, a numerical model for simulating two-dimensional, non-hydrostatic flow over irregular bottom topography has been developed. The implementation and performance of the model have been validated. The primary purpose of this model is to study flows characterized by waves or buoyancy-induced currents interacting with non-uniform topography. Applications of the model are limited to spatial domains with variability in two dimensions, which is primarily due to the high demands on the computational capacity for large scale modeling when the nonhydrostatic pressure is considered. With the developed numerical model, three numerical examples have been presented. One is the internal wave generation by the interaction of tide current with the bottom topography and its subsequent disintegration. The second is the interaction of large amplitude internal solitary wave with slope

712

topography, and the experiment setup is close to the observation in the northern continental shelf of the South China Sea, and the last one is on the stratified flow over Gaussian topography. Probably the most significant limitation of the present model is the lack of quantitative confidence in the closure scheme. However, it has been implemented in terms of spatially and temporally variable eddy viscosity and diffusivity. For DNS (i.e., Direct Numerical Simulation (DNS), at low to moderate Reynolds number), these coefficients are fixed and sufficient resolution is necessary to resolve the dissipative process. When Large Eddy Simulation (LES) is used, the values of these coefficients are determined by parameterization scheme. Other higher-order evolution equations on turbulence closure scheme have not included in the present model. Turbulence closure in stratified fluids, particularly when internal waves are considered, is an active area of current research [3,4,5,17]. However, the numerical code is based on the modular techniques of FORTRAN 95, which makes the model more flexible for further improvement.

[6]

[7]

[8]

[9]

[10]

[11]

[12]

Acknowledgement This work was supported by the PRIC Innovation Foundation. References [1] [2] [3]

[4]

[5]

FARMER D. M., FREELAND H. J. The physical oceanography of fjords [J]. Prog. Oceanogr., 1980, 12(2): 147-219. MOUM J. N., SMYTH W. D. The pressure disturbance of a nonlinear internal wave train [J]. Journal of Fluid Mechanics, 2006, 558: 153-177. BERNTSEN J., XING J. X. and GUTTOM A. Assessment of nonhydrostatic ocean models using laboratory scale problems [J]. Continental Shelf Research., 2006, 26(12): 1433-1447. SCOTTI A., BEARDSLEY R. and BUTMAN B. On the interpretation of energy and energy fluxes of nonlinear internal waves: an example from Massachusetts Bay [J]. Journal of Fluid Mechanics, 2006, 560: 103-112. LAMB K. G. Numerical experiments of internal wave generation by strong tidal flow across a finite amplitude

[13]

[14]

[15]

[16]

[17]

bank edge [J]. J. Geophys. Res., 1994, 99(C1): 843-864. HIBIYA T., OGGSAWARA M. and NIWA Y. A numerical study of the fortnightly modulation of basin-ocean water exchange across a tidal mixing zone [J]. J. Phys. Oceanogr., 1998, 28(6): 1224-1235. CUMMINS P. F. Stratified flow over topography: Time-dependent comparisons between model solutions and observations [J]. Dynam. Atmos. Oceans., 2000, 33(1): 43-72. VLASENKO V., STASHCHUK N. and HUTTER K. Baroclinic tide: Theoretical modeling and observational evidence [M]. London: Cambridge University Press, 2005. PACANOWSKI R., PHILANDER S. G. H. Parameterization of vertical mixing in numerical models of tropical oceans [J]. J. Phys. Oceanogr., 1981, 11(11): 1443-1451. STACEY M. W., ZEDEL L. J. The time-dependent hydraulic flow and dissipation over the sill of observatory inlet [J]. J. Phys. Oceanogr., 1986, 16(6): 1062-1076. NAKAMURA T., AWAJI T. and HATAYAMA T. et al. The generation of large-amplitude unsteady lee waves by subinertial K1 tidal flow: A possible vertical mixing mechanism in the Kuril Straits[J]. J. Phys. Oceanogr., 2000, 30(11): 1601-1621. CAI Shu-qun, LONG Xiao-min and HUANG Qi-zhou. Preliminary numerical study on the generation condition of induced internal solitons in the northern South China Sea [J]. Acta Oceanologica Sinica, 2003, 25(4): 119-124 (In Chinese). CAI Shu-qun, GAN Zi-jun and LONG Xiao-min. Some characteristics and evolution of the internal solitons in the Northern South China Sea [J]. Chinese Science Bulletin, 2002, 47(1): 21-26. XU Zhao-ting, YAO Feng-chao and WANG Mei. Generation of nonlinear internal waves on continental shelf [J]. Journal of Hydrodynamics, Ser. B, 2001, 13(3): 127-132. XU Zhao-ting, SHEN Guo-jing and QIAO Fang-li. Fission laws of initially interface solitary waves in two layer ocean [J]. Journal of Hydrodynamics. Ser. B, 2004, 16(5): 548-554. ORR M. H., MIGNEREY P. C. Nonlinear internal waves in the South China Sea: Observation of the conversion of depression internal waves to elevation internal waves [J]. J. Geophys. Res., 2003, 108(C3): 3064, doi:10.1029/2001JC001163. PRASAD D. Dynamics of large amplitude internal solitary wave in stratified flow over topography [D]. Ph. D. Thesis, Cambridge, USA: Massachusetts Institute of Technology, 1997.