京都大学防災研究所年報 第 56 号 B 平成 25 年 6 月
Annuals of Disas. Prev. Res. Inst., Kyoto Univ., No. 56 B, 2013

Density Current Simulation Using Two-Dimensional
High Resolution Model

M. R. SOLIMAN(1) , Satoru USHIJIMA(2) , Nobu MIYAGI(1) and Tetsuya SUMI
(1) Graduate School of Engineering, Kyoto University(2) Professor, Academic Center for Computing and Media Studies, Kyoto University
Synopsis
A two-dimensional multi-phase numerical model for incompressible, immiscible and variable density fluids has been developed based on Navier-Stokes equations. The governing
equations were discretized in a collocated grid system and solved with a finite volume
method. The present model was parallelized in a shared memory system using openMp
which reduced the computational time. The model was verified against Boussinesq and nonBoussinesq lock exchange problems. In addition, the model was applied to simulate salinity
intrusion into Ohashi River that connects Nakaumi and Shinji brackish coastal lakes. The
developed model successfully predicted the density current propagation in the experimental
cases while it still needs refinement for field study cases.

Keywords: Numerical modeling, Multi-phase flow, Density current, OpenMPtion. This stratification acts as a barrier restraining the

1. Introduction

mixing of the water column in coastal lakes. As a reGravity or density currents constitute a large class ofsult of the incomplete mixing of the water column andnatural flows that are generated and driven by the den-lack of light for the photosynthesis in deeper layers,sity difference between two or even more fluids. The

the water column may become anoxic. Assessing thedensity difference between two fluids usually arises

exchanges between water layers is important because

due to differences in temperature or salinity, but it canthey influence the trophic organization of the ecosys-

also arise due to the presence of suspended solid par-tem. Density stratification in lakes, estuaries and inticles. In water environments, saline intrusions or oil

the oceans leads to circulation and mixing of large vol-spills in the oceans, turbid water intrusions in a lakeumes of water resulting in multiphase flow. Modelingand suspended sediment plumes in a river are good ex-of these processes is of great importance for conserva-

amples. In case of shallow coastal lakes or estuaries,

tion of the ecological balance of these systems.

usually there is an interaction between sea water andSeveral laboratory experiments have been performed

the lake’s fresh water. During tide, the sea water prop-

to study density currents (e.g. Lowe et al., 2005; Kolar

agates into the lake, causing the sea water to plunge

et al., 2009). With the rapid development of numerical

underneath the free surface of a receiving water basin.

methods and advancements in computer technology,

The plunge sea water forms a density current that con-Computational Fluid Dynamics has been widely usedtinues to flow along the bottom causing water stratifica-

to study density currents. Mathematical and numeri-― 591 ―

∫

I
∫
∂ui
dΩ +
ui u · n dl =
fi dΩ
Ω ∂t
∂Ω
Ω
I
I
1
−
p · n dl +ν ∇ui · n dlρ ∂Ω∂Ωcal models when properly designed and tested against
field or laboratory data, can provide significant knowledge for density current dynamics and characteristics.
The characteristics of a gravity current head have been

(2)

studied by means of Large-Eddy Simulation (LES) andwhere t is time; ρ is the cell averaged fluid density; u is

3D Direct Numerical Simulations (DNS) of flow fronts

the velocity vector; ui is the velocity component within the lock-exchange configuration (e.g. Hartel et al.,subscript i = 1 and 2 in x- and z- directions, respec-

2000, Ooi et al., 2007). Other researchers studied thetively; x and z are the vertical orthogonal coordinates;

density current using commercial CFD codes basedn is the unit normal vector on the domain boundary ∂Ωon multiphase flow as FLUENT and FLOW-3D (e.g.

towards the outside direction; p and ν are the average

Georgoulas et al., 2010; Sangdo et al., 2012). A num-

fluid pressure and dynamic viscosity, respectively; andber of hydrodynamic models are also used to simu-

fi is the external force component due to acceleration.

late density currents in lakes, estuaries and coastal seas(e.g. Kasem and Imran, 2004; De Cesare et al., 2006;

2.2 The interface tracking method

Ushijima et al., 2007)

The motion of the interface between immiscible liq-

In this paper, a 2D vertical multiphase model has

uids of different densities and viscosities is defined bybeen developed for studying salinity propagation pro-the volume fraction function φ and the interface is con-cess in brackish coastal lakes due to tidal water flux

vected using the computed velocity field by the follow-

fluctuations. The present model describes a multi-

ing three conditions:

phase system for immiscible and incompressible flu
= 1,
fluid 1






φ=
0 < > 1, interface






= 0,
fluid 2

ids with collocated grid system which are parallelized
in a shared memory system by Open Multiprocessing(OpenMp, 2011). The present research aims to present
the validity, usefulness and applicability of multiphase
numerical approach for the simulation and study of thehydrodynamic and salinity intrusion that are usually(3)

The time dependence of the volume fraction func-

the density current propagation. The developed modeltion φ is governed by:
∫
I
∂φ
dΩ +
uφ · n dl = 0Ω ∂t∂Ωis validated against the experimental data for Boussi-

Eq.(4) states that φ moves with the fluid and is consid-

nesq and non-Boussinesq lock-exchange problems. In

ered as a flag identifying cells that contain fluid. Ac-

addition, saltwater intrusion into the tidal Ohashi Rivercording to the local value of φ, appropriate properties

is simulated and validated using field measured data.

and variables are assigned to each cell within the do-formed at estuaries and coastal lakes due to tide effectrepresenting the water surface fluctuations as well as

(4)

main.

2. Numerical Procedures

2.3 Convection diffusion equation for transported
concentration

2.1 Flow governing equations

The transported concentration is obtained by solving

The basic governing equations for flow are themass conservation equation and momentum equations,Eq.(1)-Eq.(2), which describe the flow conditions in
the vertical domain Ω as follows:
∫

∂ρ
dΩ +
Ω ∂t

∂Ω

(5)

where c is the fluid concentration, and σc is SchmidtIρu · n dl = 0

the convection-diffusion equation:
∫
I
I
∂c
ν
dΩ +
uc · n dl =∇c · n dl
∂t
σ
Ω
∂Ω c
∂Ω

number for concentration diffusion, which represents(1)

the ratio of viscosity to fluid diffusivity.

― 592 ―

two-phase immiscible fluids (liquid and gas); ν is

2.4 Computational procedure
The Numerical computation consists of four stages:

the cell average viscosity according to the existingflow computation, interface tracking, mass transport,

fluid; νl and νg are the viscosities of the two-phaseand updating fluid properties stages. The governing

immiscible fluids; ρo and co are the reference den-equations are discretized in a collocated grid systemsity and concentration values which are usuallyand solved with a finite volume method.

taken as 0.998 g/cm3 and zero, respectively, forfresh water; and β is the direct proportional coefficient of the salinity gradient between the maxi-(i) In flow computation stage; the flow equations,

mum and minimum value.

Eq.(1)-Eq.(2), are solved according to the existingdomain and boundary conditions. Estimated ve-(v) Finally, the model goes to the flow computation

locity components are obtained at the cell center

stage again to calculate the water pressure and theusing pressure at the present computational step

flow velocity field and so on.

while upwind scheme is used to calculate convective terms. The cell center velocity components

2.5 Open multiprocessing computation technique(OpenMP)are interpolated at the cell boundaries. C-HSMAC
algorithm is used to correct the pressure field so

A parallel computation technique in shared mem-that the divergence of velocity in a cell should be

ory system is implemented in the present model. Thisless than a small threshold value (Ushijima et al.,model is parallelized by Open Multiprocessing (Open2003). This algorithm for pressure computation

MP, 2011). OpenMP is an implementation of multi-

has the power to deal with large density contrast

threading which is a method of parallelizing whereby

fluids (Ushijima et al., 2007). The cell centereda master thread (a series of instructions executed con-velocity components then are corrected with the

secutively) forks a specified number of slave threads

obtained pressure values.

and a task is divided among them, Fig.1. The threads(ii) In the interface tracking stage; the volume frac-then run concurrently, with the FORTRAN runtime en-tion function φ is calculated from Eq.(4) using the

vironment allocating threads to different processors.

computed velocity field. The interface between

This process enables us to significantly reduce thedifferent fluids is achieved and hence various fluid

elapsed time of the large-scale computations.

zones can be determined.
(iii) In the mass transport stage; concentrations canbe estimated using the convection-diffusion equa-

3. Model Validation With Experimental Resultstion, Eq.(5), using the computed velocity field.
(iv) In the updating fluid properties stage; a simpli-

The present model has been tested to evaluate its ap-

fied equation of state is used to compute fluid

plicability. Therefore, it was applied to density current

density according to the volume fraction function

problems and compared with experimental results andand concentration values, given in a simple linear

field data. The simulation results for these test cases

equations:are described below.

ρ = φρl + (1 − φ)ρgν=

φρl νl + (1 − φ)ρg νg
ρ(6)

3.1 Lock-exchange problem
The lock-exchange test or dam break problem is de-(7)

fined as water of different densities are separated bya vertical barrier which is removed at time zero. Lockwith
ρl = ρo + β(c − co )(8)exchange involving fluids either with small density differences (the Boussinesq case) or with large density

where ρ is the cell average density according todifferences (the non-Boussinesq case). The removal

the existing fluid; ρl and ρg are the densities ofof the barrier results in a gravity current of the lighter― 593 ―

Fig .1 An illustration of multithreading process where the master thread forks off a number of threads which
execute blocks of code in parallel.

length (B) = 58.4 cm, height (H) = 29.5 cm, and width= 2.54 cm. The both mixing fluids were initially separated by the barrier which was 0.4 mm thick of stainless steel sheet. The both mixing fluids were fresh water and saline water with a salt concentration of 17.5
ppt. Water at room temperature was used, resulted in aheavy water density (ρ1 ) of = 1.011 g/cm3 and a light

Fig .2 Sketch of the lock-release tank.

water density (ρ2 ) of 0.998 g/cm3 , for a density ratio of
0.987, which was clearly in the Boussinesq realm.

fluid propagating at constant speed along the upper sur-Herein, we present high resolution outcomes for di-face of the channel into the heavy fluid. In the op-

rect and quantitative comparison of experimental andposite direction a gravity current of the heavier fluid

numerical results throughout the domain, not just the

propagates also at constant speed along the bottom ofwave front. Fig.3 compares the model results with thethe channel. The lock-exchange test represents mix-laboratory experiments and determines how well theing processes that occur frequently in geophysical fluid

model captured the evolution of a density field that isdynamics. For the Boussinesq case: when a fresh wa-

driven entirely by gravity currents. This figure clearly

ter river empties into a salt water estuary while for theshows a good agreement between the calculated den-non-Boussinesq case: the volcanic eruptions often takesity propagation of the wave front compared with thosethe form of gravity currents. The density within thewhich are experimentally measured. In addition, a co-flow is a result of suspended ash and hot rocks, andincidence between the front position for the estimatedis often many times larger than the surrounding air.

and the measured data was obtained as shown in Fig.4.

Fig.2 presents a sketch for the used the lock-release
cell, showing two fluids of densities ρ1 and ρ2 , eachhaving a depth of H and separated by a removable

3.3 Non-Boussinesq lock-exchange experiment

lock gate. Herein, the Boussinesq and Non-Boussinesq

Lowe et al. (2005) have carried out a series of exper-

lock-exchange tests have been used to validate the de-iments on the lock exchange involving fluids with large

veloped code.

density differences - the non-Boussinesq case. One of
the experiments is selected in this study. The experi-

3.2 Boussinesq lock-exchange experimentments were performed in a rectangular channel using

This model was applied to the Boussinesq lock-

two fluids of density ratio (γ = 0.681). The less dense

exchange experimental results obtained by Kolar et al.

fluid was freshwater (ρ2 = 0.998 g/cm3 ) and the denser(2009). The density cell had internal dimensions of

fluid was either a solution of sodium chloride (NaCl)― 594 ―

Fig .3 Comparison of experimental and predicted results for the density current propagation sequences for the Boussinesq lock-exchange case;(left) laboratory results and (right) model results.

Fig .5 Comparison of experimental and predicted results for the density current propagation sequences for the non-Boussinesq lock-exchange√
case, (where t∗ = t (g(1 − γ))/H); (left) model
results and (right) laboratory results.

or sodium iodide (NaI) (ρ1 = 1.466 g/cm3 ). The channel length was B = 182 cm, 23 cm wide and was filled
to a depth of H = 20 cm. The flow was started by avertically removing lock gate.
Fig.5 compares the model results with the laboratory experiments and determines how well the model
captured the evolution of a density field that is driven
entirely by gravity currents. This figure clearly showsan acceptable agreement between the calculated density propagation of the wave front compared with thosewhich are experimentally measured but in a different

Fig .4 Comparison of the front position estimatedfrom the laboratory experiments and the present
model for the Boussinesq lock-exchange case.

time scale. In addition, the developed model successfully predicted the bore within the heavy front. The
front velocity of the density current is almost the samefor the light case as well as the heavy case as shownin Fig.6. In the heavy front case, discrepancy may be― 595 ―

Fig .6 Comparison of the relative front position estimated from the laboratory experiments andthe present model for the non-Boussinesq lockexchange case; (a) light front and (b) heavy front

Fig .8 Initial conditions at 10:30 Am: (a) water level
fluctuation, (b) measured salinity values along
the Ohashi River, and (c) estimated salinity values along the Ohashi River.

observed because of the gate removal effect.

4. The Model Application on Ohashi River, JapanThe Ohashi River in southwest Japan is a tidal
river located between two brackish lakes, Shinji andperature on density values. The computational domain

Nakaumi, Fig.7. These brackish lakes are stably strat-

was separated vertically into two phases, air at top andified due to salinity (density) differences. Nakaumi

water at bottom with different concentrations. A com-

Lake is strongly affected by the salt water entering

putational time step of 0.1 sec. and a computationalfrom the Sea of Japan through Sakai channel while

grid of 0.25m and 100m in the vertical and the horizon-

the freshwater Hii River influenced on Lake Shinji.

tal directions, respectively, were used in the numerical

The density gradient has a significant impact on water

calculations. The computations were performed using

movement in the Ohashi River. In this study, the field

32 threads which could reduce the computational timemeasured data (September 10, 1999) for the saltwaterby 85% compared with the case of using single thread.

intrusion into the Ohashi River are used to validate the

The linear wave theory was used to calculate the in-present model (Moriwaki et al., 2003). Fig.7 shows thetrusion water flux from Lake Nakaumi using the waternine measuring stations along the Ohashi River.

level difference between the two lakes. The model re-

The salinity initial conditions of the salinity con-sults show good agreement against the measured data

centration values, computational domain, and the twofor salinity concentrations at 14:30 Pm and 16:30 Pm

lakes water levels are introduced in Fig.8. The model

as shown in Fig.9 and Fig.10. However, model refine-run starting time was 10:30 Am. The used simulation

ment is needed to avoid the stripping shape of salinity

time period was chosen to decrease the effect of tem-concentrations.

― 596 ―

Fig .7 Geographical location of the Ohashi River and sampling locations: (a) location in the Lake Shinji-Nakaumi
system and (b) nine sampling localities in September, 1999.

Fig .9 Comparison of the measured and the estimatedresults at 14:30 Pm: (a) water level fluctuation,(b) measured salinity values along the Ohashi
River, and (c) estimated salinity values along the
Ohashi River.

Fig .10 Comparison of the measured and the estimatedresults at 16:30 Pm: (a) water level fluctuation,(b) measured salinity values along the Ohashi
River, and (c) estimated salinity values along
the Ohashi River.

― 597 ―

tion of baroclinic shallow water models: The lock-5. Conclusion

exchange problem, Journal of Ocean Modelling,
A two-dimensional multiphase numerical model forincompressible, immiscible and variable density flu-

Vol.28, pp.137-152.
Lowe, J., James, Rottman, W. and Linden, F. (2005):

ids has been developed to study salinity propagation

The non-Boussinesq lock-exchange problem: Partof brackish coastal lakes due to tidal water flux. The

1. Theory and experiments, Journal of Fluidgoverning equations are discretized in a collocatedgrid system and solved with a finite volume method.
Using openMp parallel programing approach reason-

Mech., Vol.537, pp.101-124.
Moriwaki, S., Fujii, T. and Fukui, K. (2003): Intrusion of high salinity water mass into Ohashi

ably decreased the computational time. The developed

River, San’in District of Japan, LAGUNA, Vol.10,

model successfully predicted the density current propagation in the lock-exchange problems for Boussinesq
and non-Boussinesq experimental cases. In addition,

pp.35-45.
Ooi, S.K., Constantinescu, G. and Weber, L.J. (2007):
2D Large-Eddy Simulation of lock-exchange

saltwater intrusion into Ohashi River was simulated

gravity current flows at high grashof numbers.

and validated using field measured data. The predicted

Journal of Hydraulic Engineering, Vol.133, No.9,

model still needs refinement to avoid concentration oscillations.

pp. 1037-1047, ISSN 0733-9429.
OpenMP (2011): OpenMP applicaton program interface, http://www.openmp.org.

References

Sangdo, A., Pierre, Y., Julien, Subhas, K. and Ve-

De Cesare G, Boillat G, Schleiss A (2006) Circu-

nayagamoorthy (2012): Numerical simulation oflation in stratified lake due to flood-induced tur-particle-driven gravity currents, Environ Fluid

bidity currents, Jour. Environ Engng., ASCE, Vol.

Mech, DOI: 10.1007/s10652-012-9251-6.
Ushijima, S., Okuyama, Y., Takemura, M. and Nezu,

132, No. 11, pp. 15081517.
Georgoulas, A., Angelidis, P., Kopasakis, K. andI. (2003): Parallel computational method for pres-

Kotsovinos, N. (2010): 3D Modelling of Turbid-

sure field in incompressible flows on 3d curvilin-

ity Currents, Environ Fluid Mech, Vol.10, pp.603-

ear coordinates, Annual Journal of Hydraulic Engineering, Vol.47, pp.385-390.

635.
Hartel, C., Meiburg, E. and Necker, F. (2000): Anal-

Ushijima, S., Makino, O. and Nezu, I. (2007):

ysis and Direct Numerical Simulation of the Flow

Three dimensional computational method forat a Gravity-Current Head: Part 1. Flow Topology

wave flows running up in a river and their flooded

and Front Speed for Slip and No-Slip Boundaries.

flows, Journal of Hydroscience and Hydraulic En-

Journal of Fluid Mechanics, Vol.418, (September

gineering, Vol.25, pp.31-39.
Ushijima, S., Yamada, S., Enkai, M., and Nezu, I.

2000), pp.189212, ISSN 1469-7645.
Kassem, A., and Imran, J. (2004), Three-dimensional(2007): Transportation of solid bodies in wave

modeling of density current. II. Flow in sinuous

flows predicted by 3D multiphase model: Pro-

confined and unconfined channels, Jour. Hydraul.

ceedings of 32nd IAHR Congress, Theme D-2b,

Res., Vol. 4 No. 2 (6), pp. 591602.

CDROM.

Kolar, R., Kibbey, T., Szpilka, C.M., Dresback, K.,
Tromble, E., Toohey, I.P., Hoggan, J. and Atkinson, J. (2009): Process-oriented tests for valida-― 598 ―

(Received June 11, 2013)