The spreading and mixing of dense gas clouds in still air

of

dense gas clouds in still air

A. P. van Ulden

1. De meetresultaten van de Thorney Island experimenten verdienen een zorg­ vuldige heranalyse.

(dit proefschrift, III.5.3)

2. Bij windstil weer mengt een zware gaswolk zich sneller met lucht, dan een gaswolk met dezelfde dichtheid als lucht.

3. Door het turbulentiemodel uit dit proefschrift aan te passen aan de stromingskarakteristieken van de stabiele atmosferische grenslaag, kan betrekkelijk eenvoudig een geostrofische weerstandswet worden afgeleid. (dit proefschrift, III.4; Nleuwstadt, 1985: A model for the stationary stable boundary layer)

4. Tussen de Lorentzkracht en de Magnuskracht bestaat een grote analogie.

5. Emissie en depositie van ammoniak behoren in een luchtverontreinigings­ model niet gescheiden te worden behandeld.

(Asman, 1987; proefschrift L.U. Wageningen)

6. Windconvergentie bevordert radioactieve neerslag bij kernrampen op tweeërlei wijze.

(Cats et al., 1987; K.N.M.I. wetenschappelijk rapport 87-1)

7. In het operationele weermodel van het Europese centrum voor middellange termijn weersverwachtingen wordt het mengend vermogen van thermiek ontoereikend behandeld.

8. Windsnelheidswaarnemingen boven land kunnen met locale terreincorrecties geschikt worden gemaakt als invoerparameter voor numerieke weermodellen. (Wieringa, 1986; Quart. J. Roy. Met. S o c , 112)

9. Natuurkundige wetten beperken de schaal en de diepte van atmosferische depressies. Economen zouden hier lering uit kunnen trekken.

10. Het storten van alle verontreinigde grond en afval op één plaats in

Nederland maakt het probleem van de bodemverontreiniging beter beheers­ baar en vergroot op termijn de skimogelijkheden in ons land.

11. Door korte termijn programmering en financiering van het wetenschappelijk onderzoek bevordert de overheid het ontstaan van een play-back-show cultuur bij onderzoekers.

T l J6*l ^

(l #w

6 W

THE SPREADING AND MIXING OF

DENSE GAS CLOUDS IN STILL AIR

Proefschrift

Ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus, prof. dr. J.M. Dirken,

in het openbaar te verdedigen ten overstaan van een commissie

door het College van Dekanen daartoe aangewezen, op donderdag 21 januari 1988 te 1^.00 uur

door

Adrianus Petrus van Uiden geboren te Leiden Theoretisch Natuurkundige

(TR diss

1605

Dit proefschrift is goedgekeurd door de promotoren:

Prof. Dr. Ir. G. Ooms Prof. Dr. J. Wieringa

This thesis is also published as Scientific Report W.R. 87-12 of the Royal Netherlands Meteorological Institute

CONTENTS

page

Voorwoord

Samenvatting

Summary ₁₀

Introduction, and problem definition ₁₃

Chapter I

1.1 1.2

Introduction to the modelling of quasi-uncoupled dense clouds

Basic principles Simple scaling "laws"

19 19 22

Chapter II : A dynamical integral model for two-dimensional gravity currents

11.1 Introduction

11.2 History and problem definition

11.3 A new bulk model for fixed volume releases 11.4 Some general model results

11.5 The momentum budget for great times 11.6 Summary and conclusions

29 29 29 32 43 46 47 Chapter III : III.1 111.2 111.3 111.4 111.5 111.6 111.7

The spreading and mixing of a dense cloud in still air

Introduction

The general model structure The momentum-integral equation The energy budget

Analysis and simulation of still-air experiments Model characteristics and sensitivity analysis Conclusions 49 49 51 54 62 68 78 94

Chapter IV : Reflections on the construction of a comprehensive

model for dense cloud dispersion 97

IV.1 Introduction 97 IV.2 Basic processes in passive cloud dispersion 97

IV.3 Modeling the intermediate phases 100

IV.U Concluding remarks 101

References 103

Curriculum Vitae 108

Chapter II is based on a previous publication:

Van Ulden, A.P., 1984: A new bulk model for dense gas dispersion: two-dimensional spread in still air. In: Atm. Dispersion of Heavy Gases and Small Particles (G. Ooms and H. Tennekes eds., Springer Verlag) 419-440.

#*

Chapter III is a synthesis and extension of the following papers: Van Ulden, A.P., 1987 A: The spreading and mixing of a dense cloud in

still air. Proceedings of I.M.A. conference on Stably Stratified Flow and Dense Gas Dispersion, Chester, April, 1986. Oxford University Press, J. Puttock, editor, in press.

Van Ulden, A.P., 1987 B: The heavy gas mixing process in still air at Thorney Island and in the laboratory. J. of Hazardous Materials, Vol. 16, pp 411-426.

Van Ulden, A.P., 1986: The relative importance of turbulence generated by a spreading dense cloud and atmospheric turbulence. Paper presented at the Second Symposium on Heavy Gas Dispersion Trials at Thorney Island, Sheffield, September, 1986.

Van Ulden, A.P., 1987 C: Some aspects of the dispersion of dense puffs. Paper presented at the Third international Symposium on Stratified Flows, Pasadena, February, 1987.

VOORWOORD

In dit proefschrift wordt onderzoek beschreven dat verband houdt met de preventie van rampen door gevaarlijke stoffen. Het onderzoek vindt zijn oor­ sprong in de interesse van het Directoraat Generaal van de Arbeid voor het zo veilig mogelijk hanteren van gevaarlijke stoffen. In het begin van de jaren zeventig werd onderzoek gestart naar de problematiek van chloor. Chloor is een giftig gas; het werd in de eerste wereldoorlog gebruikt als strijdgas. Het wordt thans in grote hoeveelheden geproduceerd als basisstof voor de chemische procesindustrie. Chloor is veel zwaarder dan lucht. Daarom bestond de vrees, dat bij een ongeval met chloor, dit zware gas als een deken het aardoppervlak zou bedekken en zich weinig met lucht zou mengen. Dit zou aanzienlijke

risico's voor de bevolking met zich mee kunnen brengen.

Over deze problematiek werd het KNMI in'1971 benaderd door de Arbeids­ inspectie. Als jong onderzoeker werd ik aan dit probleem gezet. Het bleek een bijzonder probleem te zijn, zowel vanwege de grote maatschappelijke relevan­

tie, als vanwege het vrijwel ontbreken van kwantitatieve beschrijvingen in de toenmalige literatuur. In de beginfase van mijn onderzoek heb ik veel steun ontvangen van Ger Abraham (Waterloopkundig Laboratorium), die mij wegwijs heeft gemaakt in de problematiek van dichtheidsstromingen. Ook heb ik enkele verkennende experimenten kunnen doen bij het Waterloopkundig Laboratorium. Op basis van deze experimenten, en op basis van in de literatuur aanwezige infor­ matie over de spreiding van olie op water, heb ik een eenvoudig model ontwik­

keld voor de spreiding van een zware gaswolk op een horizontaal oppervlak. Ondermeer met behulp van dit model werden grootschalige experimenten ontworpen en uitgevoerd door een twintigtal bedrijven en overheidsinstellingen in

Nederland. Deze experimenten werden uitgevoerd op de Maasvlakte bij Rozenburg in 1973 en 197^ onder de bezielende leiding van de heer E.H. Siccama (Arbeids­ inspectie). Ik heb voor het KNMI aan deze experimenten deelgenomen, daarbij voortreffelijk ondersteund door Wim Schipper en door de Instrumentele Afdeling van het KNMI. Voorts leverde de weerkamer in De Bilt speciale weersverwachtin­ gen, die essentieel waren voor het succes en de veiligheid van de experimenten.

De Nederlandse experimenten leverden zowel direct bruikbare informatie op, als ook inspiratie voor verder onderzoek. Eind jaren zeventig kwam dit onderzoek internationaal sterk op gang. Dit leidde tot een aantal nieuwe expe­ rimenten in het buitenland en tot een toename van de internationale samenwer­ king op het gebied van modelontwikkeling. Deze internationale ontwikkelingen

vormen de voedingsbodem voor het recente onderzoek dat beschreven is in dit proefschrift.

Zo'n 15 jaar zijn voorbijgegaan sinds de start van het onderzoek naar de verspreiding van zware gassen. Nil vordert het internationale onderzoek met rasse schreden. Bevredigende en praktisch bruikbare modellen liggen thans in het verschiet. De boodschap is duidelijk: het bereiken van tastbare resultaten met behulp van resultaat gericht onderzoek kan een aanzienlijke tijd vergen. . Onderzoeken is vooruitzien.

Ik ben het KNMI erkentelijk, dat het mij de ruimte gaf: 15 jaar geleden voor het opzetten van dit onderzoek en de laatste jaren voor het realiseren van dit werkstuk. Voorts gaat mijn dank uit naar de vele personen in binnen-en buitbinnen-enland die ebinnen-en rol hebbbinnen-en gespeeld bij het tot stand kombinnen-en van dit proefschrift. In het bijzonder ben ik dank verschuldigd aan Bronno de Haan, Theo van Stijn en Roland Stull voor het slechten van numeriek-wiskundige

barrières en aan Günther Können voor zijn spirituele inbreng. Tenslotte spreek ik mijn waardering uit over de grote inzet van Marleen Kaltofen en van de tekenkamer en drukkerij van het KNMI bij het feitelijk realiseren van dit proefschrift.

SAMENVATTING

In onze industriële samenleving worden vele gassen gebruikt en geprodu­ ceerd, die brandbaar, explosief, radioactief of giftig zijn. Dit maakt onge­ vallen mogelijk waarbij grote hoeveelheden van een gevaarlijk gas in de atmos­ feer kunnen vrijkomen. Zulke ongevallen hebben plaats gevonden en zullen met grote waarschijnlijkheid ook in de toekomst plaats vinden. Een speciale klasse van ongevallen is die waarbij in korte tijd een wolk zwaar gas vrijkomt. Zo'n zware gaswolk stroomt uit over het aardoppervlak onder invloed van de zwaarte­ kracht. Dit kan resulteren in hoge gasconcentraties nabij de grond over grote oppervlakken, met alle risico's van dien. In het geval van een zeer grote en dichte gaswolk en lage windsnelheden zal deze uitstroming weinig worden beïn­ vloed door de atmosferische omstandigheden: de gaswolk bepaalt zijn eigen ont­ wikkeling en vormt hierbij een zogenaamd quasi-ontkoppeld systeem.

Dit proefschrift is gericht op de beschrijving van quasl-ontkoppelde gaswolken. De opbouw van het proefschrift is als volgt.

In een korte inleiding worden de belangrijkste vormingsmechanismen van zware gaswolken gegeven en het verloop van de verspreiding geschetst. Dit wordt geïllustreerd met een aantal foto's.

Hoofdstuk I beschrijft enkele elementaire modelleringsprincipes. Er wordt een aantal eenvoudige "schalingswetten" afgeleid, die ruwe schattingen mogelijk maken van de afmetingen van en de concentraties in zware gaswolken. In hoofdstuk II wordt de problematiek van twee-dimensionale dichtheids­ stromingen geanalyseerd. 'Aandacht wordt besteed aan de structuur van de voorste begrenzing van de zware vloeistof (of het zware gas). Zo'n voorste begrenzing wordt gekenmerkt door een abrupte sprong in de dichtheid en wordt veelal front genoemd. Twee-dimensionale fronten zijn goed bestudeerd in het laboratorium. Experimentele informatie uit de literatuur wordt in dit hoofd­ stuk gebruikt voor het afleiden van schalingswetten voor de frontstructuur. Voorts wordt in dit hoofdstuk een model afgeleid voor de dynamica van de zware vloeistof achter het front. Hierbij wordt gebruik gemaakt van ruimtelijk geïn­ tegreerde bewegingsvergelijkingen. Een dergelijk model wordt dynamisch

integraal-model genoemd.

Hoofdstuk III vormt de kern van dit proefschrift. Hier wordt een drie­ dimensionaal dynamisch integraal-model voor een zware gaswolk afgeleid. Het

inzakken, uitspreiden en opmengen van de wolk wordt beschreven door ruimtelijk geïntegreerde bewegingsvergelijkingen. Hierbij wordt rekening gehouden met

horizontale en verticale versnellingen in en rond de wolk. Voor de cirkelvor­ mige rand van de wolk wordt de frontstructuur van hoofdstuk II gebruikt. Het opmengen van de wolk wordt beschreven middels een ruimtelijk geïntegreerde vergelijking voor de turbulente energie. Met deze vergelijkingen worden de straal van de wolk en de gasconcentraties berekend als functie van de tijd. Deze worden vergeleken met experimentele gegevens uit de literatuur. Het blijkt dat de straal van de wolk en de concentratieprofielen goed worden beschreven. Tevens blijkt uit de evaluatie van de meetgegevens dat concentra­ ties nabij de grond vele malen hoger zijn dan tot dusverre werd aangenomen (III.5.3).

Het hoofdstuk wordt afgesloten met een uitvoerige evaluatie en gevoelig­ heidsanalyse van het model en met een beschouwing over de toepasbaarheid van simpele schalingswetten.

Ter afronding van het proefschrift wordt in hoofdstuk IV gesproken over het opzetten van een model, dat naast de quasi-ontkoppelde toestand ook andere fasen in het verspreidingsproces kan beschrijven. In zo'n model begint de wolk in de quasi-ontkoppelde fase, spreidt zich uit en evolueert uiteindelijk naar de zogenaamde passieve eindfase. In het hoofdstuk wordt in het kort beschreven welke stadia de wolk hierbij doorloopt en wat de dominerende processen zijn, welke beschrijvingen hiervoor reeds voorhanden zijn en welke problemen nog moeten worden opgelost. De conclusie is dat het dynamisch integraal-model uit hoofdstuk III een geschikt uitgangspunt is voor het construeren van een derge­ lijk algemeen model.

SUMMARY

The present work was originally started in the early seventies. In the Netherlands the authorities and companies questioned the safety of production, transport and processing of large quantities of chlorine. Chlorine is a toxic and dense gas. A release of a large amount of chlorine is a considerable hazard to the environment. However, the details of the spreading mechanism of high density gases in the atmosphere were greatly unknown. In 1973 experiments were carried out in the Netherlands to investigate these problems. Pictures, of these experiments, shown in the introductory section of this thesis visual­ ize the spreading of a dense gas cloud. These pictures show that a dense cloud spreads over the ground and takes the form of a large shallow pancake. The dominating mechanism is the action of gravity. Indeed, in the life cycle of a dense cloud there is a period in which gravity-induced spreading and mixing almost fully determines the dispersion of the cloud. So the cloud behaves as if it is uncoupled from the atmospheric flow around it. Therefore we call this phase the uncoupled phase. This work focusses on the modeling of quasi-uncoupled dense clouds. The modeling is described in three steps, subsequently given in the chapters I, II and III.

In chapter I the basic modeling principles are described and some scaling laws are derived. These allow a simple first-order estimate of the size and mean concentration of the cloud as a function of time.

In chapter II a detour is made to investigate two-dimensional density currents. The frontal structure of the leading edge of such currents has been studied relatively well in the laboratory. Laboratory data are used to derive a parameterisation of this frontal structure. The bulk dynamics of

two-dimensional density currents is studied by means of a special integral-form of the momentum equation that includes the effects of entrainment and allows the existence of large density differences between the dense flow and the ambient fluid.

In chapter III a dynamical integral-model for axisymmetric clouds is presented. From the fundamental equations of motions, integral-equations are derived for the spreading, slumping and mixing of the cloud. These equations account for radial and vertical accelerations in and around the cloud and for the effect of large density differences between the cloud and the environment. Turbulent mixing is described with an entrainment equation, which uses the turbulent kinetic energy of the cloud. This turbulent energy is described with

a newly developed time-dependent turbulent-energy equation.

The model is compared with observations on axisymmetric dense clouds taken from the literature. First radial gravity spreading is considered

(III.5.2). Here the agreement between model simulations and experimental data is quite satisfactory. The present results on radial spreading are consistent with the results of previous studies. However in this study no ad-hoc fitting of model coefficients was needed to achieve this.

Next experimental concentration data are analysed. The present analysis shows that near-surface concentrations are likely to be much higher than indicated by previous studies of the same experimental data. The reason for this is that usually a Gaussian or uniform profile is assumed for the vertical distribution of area-averaged cloud concentrations. The present analysis shows that such profiles are poor approximations to the observed profiles, especial­ ly close to the ground where strong vertical gradients are observed. In this study a new ad-hoc similarity profile is proposed. With this profile our model gives a satisfactory simulation of observed concentrations (III.5.3).

The remainder of chapter III is devoted to a sensitivity analysis and to a description of model characteristics.

In chapter IV we put the present work in a wider context. We give our views on the possibilities to construct a satisfactory comprehensive dense gas model that includes the transport and dispersion by wind and atmospheric

turbulence. Our conclusions are that such possibilities exist and that our dynamic model for still air is a good starting point.

INTRODUCTION AND PROBLEM DEFINITION

In modern society many gases are used and produced which are flammable, explosive, radioactive or toxic. Releases of such gases have occurred in the past and are likely to occur in the future. A brief review of recent accidents has been given by Graziani et al. (1987). It appears that accidental releases of hazardous gases are a grave problem. This is especially the case when the release leads to the quasi-instantaneous formation of a large cloud. In a significant fraction of (potential) accidents, such a cloud will have a density that is larger than the atmospheric air density. Such a high density may result from the high molecular weight of the gas that is released, or from the low temperature at which the gas is released. The latter will occur when the gas has a low boiling point and is stored in liquefied form. Three major release mechanisms may result in the quasi-instantaneous formation of a dense cloud.

Type 1

The gas is stored in liquefied form under pressure at ambient tempera­ ture. After a failure of the storage tank the pressure drops suddenly and a significant fraction of the liquefied gas evaporates spontaneously. In this case the heat source for the evaporation is the heat stored in the liquefied gas. Spontaneous evaporation leads to a rapidly expanding hemispheric cloud of a dense gas. The rapid expansion causes vigorous initial mixing with the

ambient air.

Type 2

The gas i s stored in liquefied form in a container t h a t i s cooled down to the b o i l i n g temperature of the gas at ambient p r e s s u r e . In t h i s c a s e no heat i s a v a i l a b l e in the liquid for spontaneous evaporation. After t h e f a i l u r e of the tank the l i q u i d w i l l spread over the ground or over a water s u r f a c e . Especially in the l a t t e r case the heat t r a n s f e r from the surface to t h e l i q u e f i e d gas w i l l lead t o rapid formation of a dense gas cloud.

Type 3

A dense gas is stored in gaseous form and is released as such (without an intermediate boiling process). In the special case in which the gas is stored at ambient pressure, expansion and enhanced mixing will not occur during the release.

These three types of release mechanisms have been simulated in field and laboratory experiments. Examples are:

type 1: Freon experiment at the "Maasvlakte" (Van Uiden, 1971*).

type 2: Low wind cases of LNG, LPG experiments at "Maplin Sands" (e.g. Colen-brander and Puttock, 1984).

type 3: Freon experiments at Porton (Picknett, 1981) and "Thorney Island" (e.g. McQuaid, 1984).

Laboratory experiments with Freon (Havens and Spicer, 1984, 1985). These and other experiments have supported and stimulated significant modelling efforts during the last 15 years.

The modelling of the behaviour of large dense clouds is by no means a simple task. Apart from the usual problems which arise, when one wants to model the dispersion of passive trace contaminants, dense clouds create a few problems of their own. The main reason for this is that the force of gravity works on dense clouds. This leads to gravity-driven flow phenomena with their own specific characteristics. For a review of this subject we refer to Simpson

(1982).

A suitable introduction to the problems related to dense cloud modelling is a time series of pictures of a dense cloud during and after its release. For this purpose we have selected a number of pictures taken during the first large scale field experiments on dense cloud dispersion, that have been

carried out in the Netherlands in 1973. A brief survey of these experiments is given in the report "Experiments with Chlorine" (Directorate-General of

Labour, 1975). The dense gas problems are discussed in Van Ulden (1974). The aim of these experiments was to acquire a better insight into the behaviour and control of large accidental chlorine spills. To avoid the hazards of a large chlorine release, a potential accident was simulated using freon-12 (CC12F2). This is an invisible, non-toxic and rather inert gas, with a

molecular weight of 121 and a boiling point of 243 K.

A quasi-instantaneous source was obtained as follows. A vessel with 1000 kg of liquid freon at ambient temperature was placed over an open tank filled with hot water. At the bottom of the vessel was an outlet which was closed by a rupture disc. The disc was forced to break by an increase of the pressure on the vessel to about 7 bar. This led to a rapid downward outflow of

the liquid freon into the hot water. The water provided sufficient heat for a complete evaporation of the freon in about 5 seconds. Moreover it supplied the cloud with water vapour, the condensation of which made the cloud visible. After its expansion the cloud had a temperature nearly equal to the air temperature. Its density was about 20$ higher than that of air.

The formation of the cloud and its subsequent behaviour is illustrated with 6 characteristic pictures on the next pages. The first picture gives the cloud after 1 s, the second picture after 5 s when the cloud formation is completed. We see that the expansion is more or less hemispheric. The cloud edge is very diffuse. This is due to the vigorous turbulent mixing owing to the rapid expansion. After its formation the cloud starts to slump. The

pictures 3~5 show the cloud after 8 s, 15 s and 35 s respectively. We see that the cloud acquires a very shallow and wide form. In other pictures, not shown here, it is observed that this form is more or less cylindrical. After 35 s the diameter of this cylinder is about 100 m. Its depth is then about 0.5 m. It is clear that a dense cloud behaves differently from a puff of passive contaminant. In passive puff dispersion the vertical and horizontal scales are of the same order.

The very pronounced slumping, which is observed for a dense cloud, is caused by the combination of the high cloud density and gravity. This makes the cloud negatively buoyant. Therefore dense gases are often called heavy gases. The negative buoyancy results in an increased static pressure in the cloud, which leads to a radially outward pressure force. This causes the outward slumping motion. The process resembles the spreading of a liquid over a solid surface.

The slumping process does not continue indefinitely. In the absence of wind and atmospheric turbulence, the slumping will stop when the cloud height reduces to the hieght of the roughness elements and irregularities of the surface. When there is wind, and turbulence, sooner or later atmospheric

transport and diffusion will carry the cloud away and dilute it to smaller and smaller concentrations until it behaves as a passive contaminant. The start of the latter processes is shown in picture 6, which was taken after 70 s. Here we see that only remnants of the visible cloud still linger at the surface. The bulk of the cloud has become invisible due to mixing with dry air and evaporation of the liquid water.

CO co IT) II CO 00 ii CM 0) - P Ü • H PM - P Ü PH +> OH

P i c t u r e k

t=35S

P i c t u r e 5

— J L ^ *

— Cloud diameter ca100m

t=70S

These and other observations suggest that the life cycle of the spreading and mixing of a dense cloud can be split up in five phases.

I The formation phase

This phase comprises the release of the (liquefied) gas and, if appropri­ ate, its following evaporation and expansion. The three major types of formation mechanisms have been given above (picture 1).

II The quasi-uncoupled slumping phase

This phase occurs when the effects of turbulence and momentum of the ambient flow can be neglected, as well as the influence of surface fric­ tion. In this phase the dynamics and mixing of the cloud are exclusively driven by the negative buoyancy of the cloud. A slumping and radially spreading cloud is observed with a more or less pronounced raised edge, sharp radial boundaries and a diffuse top. Significant mixing occurs, resulting from the turbulence which is generated by gravity spreading (pictures 2-5).

III The coupled spreading phase

In this phase gravity spreading still dominates over horizontal atmos­ pheric diffusion. However, the cloud no longer slumps because atmospheric turbulence leads to vertical mixing at a similar or greater rate than the gravity induced slumping motion. In addition, advection by the mean wind occurs and surface friction may affect the gravity spreading (picture 6).

-IV The mixed phase

In this phase gravity spreading and horizontal diffusion are of comparable importance. The sharp-edged cloud transforms into a diffuse cloud.

Vertical diffusion is still reduced by cloud density.

V The passive phase

In this phase cloud density effects can be neglected. Dispersion by atmos­ pheric turbulence and advection by the mean wind are similar to those for a passive puff.

The present work focusses on the second (quasi-uncoupled) phase and in particular with the special case of an isothermal dense cloud, i.e. a cloud with a uniform temperature equal to the air temperature.

Chapter I

INTRODUCTION TO THE MODELLING OF QUASI-UNCOUPLED DENSE CLOUDS

1.1 Basic principles

It appears from observations that a dense cloud in the quasi-uncoupled phase has a cylindrical shape. The radial cloud edge is quite distinct (figure 1.1). This calls for a description of the cloud area involving the distance R of the radial cloud edge from the cloud center. The cloud area is given by

A = TTR2 . (1.1)

On the other hand the cloud top i s quite diffuse and i r r e g u l a r (figure 1.1). Clearly the v e r t i c a l d i s t r i b u t i o n of dense material in the c y l i n d r i c a l cloud i s non-uniform. This c a l l s for an i n t e g r a l d e f i n i t i o n of cloud h e i g h t . The mean height <z> of dense material in the cloud i s by d e f i n i t i o n

<z> = fn ƒ" z Ap(r,z) 2irr dr dz / /R ƒ" Ap(r,z) 2irr dr dz (1.2)

o o o o

Here z i s the height above t h e ground, r the r a d i a l d i s t a n c e from t h e cloud center and Ap(z,r) the difference between the l o c a l cloud density p ( z , r ) and

the ambient a i r d e n s i t y p . The mean height <z> can be used t o define a mean

3

cloud depth H. Let us consider the special case in which initially the cloud is a cylinder with radius RQ, height HQ and a uniform density difference Ap .

For this special case it can be simply shown that the relation between the cloud top HQ and the mean height <z>0 is

H = 2<z> . (1.3) o o

This suggests that a suitable general definition for the cloud depth is

H ■ 2<z> . (1.4)

Adopting this definition we may also define a cloud volume by

H

77777^77777777777777777777777777777777777777777777"

4 R M R ►

figure 1.1: Vertical Cross-section through the cloud center showing the sharp

cloud edge and the diffuse cloud t o p .

a mean cloud density difference by

R m

Ap = ƒ ƒ A p ( r , z ) 2irr dr dz / V ( 1 . 6 )

and a mean cloud density by

^p" = pa + Ap" . (1-7)

These definitions have the important and convenient property that they allow a simple and exact representation of the potential energy of the cloud. This potential energy PE equals

P = fR r g Ap(r,z) z 2irr dr dz , (1.8)

E 0 0

where g is the acceleration by gravity. The integral in (1.8) can be simply evaluated by using (1.2), (1.1), (1.5) and (1.6). The result is

-PE-.-fc g Ap" V H . (1*9)

In this equation g Ap V is the total negative buoyancy of the cloud. Thus the potential energy equals the negative buoyancy multiplied by the height ^H of the center of mass of the cloud. This result can be simplified further for the special case of an isothermal cloud. We use the fact that the total mass of material should be conserved. This mass equals p V . Moreover the total mass of the cloud increases due to turbulent entrainment of air. For an isothermal cloud the mass increase at a given ir

cloud mass at a given instant equals

cloud the mass increase at a given instant equals p (V-V ) . Thus the total 3. O

p V = p V + p (V-V ) . (1.10)

From t h i s equation i t follows d i r e c t l y that

Ap V = AJT V . (1.11) o o

Thus for an isothermal cloud the difference between the cloud mass and the mass of an amount of air with the same volume is a conserved quantity. There­ fore the negative buoyancy g Ap V is also a conserved property for an isother­ mal dense cloud. Using this we may write for the potential energy

PE = ^ g ^ V o H ' ( I J 2 )

So isothermal clouds have the convenient property that their potential energy is fully specified by the initial conditions and by the actual cloud height H.

We will use this property now to illustrate the cloud dynamics in terms of its energy budget. This may help to understand the basic principles of the spreading and mixing processes.

As initial condition we take a still cloud with a potential energy

P E o = * g ^ V o H o - ( I'1 3 )

The negative buoyancy, which is the basis of this potential energy, also creates an increased static pressure in the cloud. This increased static pressure pushes the cloud radially outward: it starts a radial acceleration, that is accompanied by a downward acceleration for continuity reasons. In terms of energy this

corresponds with a transformation of potential energy into kinetic energy (Kg). Indeed slumping implies a loss of height and a loss of potential energy.

In the slumping process the outward moving cloud edge feels the resis­ tance of the ambient air. Continually, new air is accelerated by the cloud edge. Here strong shears in the flow are present that lead to the production of turbulent energy (TE) at the cost of kinetic energy.

The turbulence thus created leads to two new processes. The first process is turbulent mixing at the top of the cloud. Ambient air is entrained into the cloud volume, which corresponds with a relative thickening of the cloud. This is equivalent with a transformation of turbulent energy back into potential energy. This process we call buoyant destruction. The second process is viscous dissipation. In a turbulent energy cascade energy is transferred to smaller and smaller scales. At the smallest scale viscous friction destroys turbulent

energy, while internal heat (IE) is produced. The flows of energy are summar­

G production of Kg S production of Tg

B buoyant destruction of Tg D viscous dissipation of Tg

In the life cycle of the quasi-uncoupled phase the potential energy decreases monotonously, while the internal heat increases. The kinetic energy starts at a zero value, increases when the cloud accelerates from rest, and decreases later on. The turbulent energy lags behind the kinetic energy, but similarly goes through a maximum during the quasi-uncoupled phase.

It will be clear that the gravity spreading of a dense cloud is a very unsteady process. A proper description of this process requires the solution of time-dependent equations for momentum and energy. This is done in chapter II for 2-dimensional gravity currents and in chapter III for axisymmetric clouds. Before we start with the detailed modeling problems, we first give a descrip­ tion of some simple scaling "laws". This may help the reader to develop some quantitative understanding of the spreading and mixing of dense clouds.

1.2 Simple scaling "laws

The source for all dense gas motions is the static pressure surplus in the cloud due to its negative buoyancy. Thus we start our scaling analysis with an evaluation of this pressure surplus. Inside the cloud the static pressure Pc

varies with height as

8P /9z' = -gp , (1.14)

where z* is the local height and p is the local cloud density. To find Pc at a

given position (r,z), we integrate (1.14) with respect to z' from z to infinity

P (r,») - P (r,z) = - ƒ gp(r,z')dz'

c c z (1.15)

Outside the cloud the static pressure Pa at the height z is given by

P (») - P (z) = - ƒ gpo dz' .

3, 3 Z a.

At an infinite height the pressure is no longer affected by the presence of the cloud, thus

P (r,-) = P (-) . (1.17)

C a

The static pressure surplus at a given position in the cloud is then given by

nQ(r,z) ■ P (r,z) - P (z) = f g Ap (r,z») dz' . (1.18)

a c ci z

A mean value n„ of the static pressure surplus may be defined in the same manner as we have defined Ap in (1.6). It reads:

no = J * ^ nQ(r,z) 2irr dr dz / V . (1.19)

Substituting (1.18) for II- we o b t a i n :

'Ü = r r g A p ( r . z ' ) 2irr dr dz dz / V . (1.20)

s o z

Using (I.2)-(I.6) it can be shown that (1.20) is identical to

1 = \ g Ap H . (1.21)

s

This result is quite general in the sense that no assumptions have been made regarding the density distribution in the cloud. However, when we want to arrive at simple scaling laws, we have to make such assumptions.

A scaling law for radial gravity spreading can be derived as follows. The radial edge of the cloud has the appearance of a sharp front (see figure 1.1). At this front there is a jump in the static pressure. We approximate this

pressure jump by its mean cloud value (1.21). The front moves outward with a velocity

U ■ dR/dt . (1.22)

The moving front meets c o n t i n u a l l y undisturbed ambient a i r which i s a c c e l e r a t e d outward and upward. This c r e a t e s a dynamic pressure n outside the cloud, which scales as

where c is an empirical drag coefficient. Assuming a balance between the mean static pressure inside the cloud and the dynamic pressure outside the cloud (i.e. nc = n ) we simply find from (1.21) and (1.23) that

o ex

Uf = k (g Ap H / pa) * (1.24)

where

k ■ 1 / /c (1.25)

is a densimetric Froude number. We assume that cloud shape and flow configur­ ation only depend on H, R and Uf. Then k is a constant. Experimental data show

that the simple scaling law (1.24) gives a fair description of the radial spreading of a dense cloud and that k is close to unity, except initially when the cloud is still accelerating from rest (see chapters II and III). From this result we may draw a first and quite important conclusion, that is that the radial spreading velocity does not depend on the density surplus Ap alone, but on the product of Ap and H! Thus significant gravity spreading may occur for small values of Ap provided the cloud depth is great enough. Another point, to be noted is, that it is the ratio Ap/p that matters and not Ap/p. The latter

3.

r a t i o has been used by several authors ( e . g . Van Ulden, 197*0. A numerical example i l l u s t r a t e s the kind of spreading v e l o c i t i e s which we may expect. Taking Ap/p = 0.1 and H = 10 m we find t h a t Uf = 3 m s "1. For the same cloud

depth and Ap/p = 1 we find Uf = 10 m s- 1.

3.

From (1.24) we may derive an expression for the cloud radius as a function of time. Using (1.5), (1.11) and (1.22), we may rewrite (1.24) as

y

dRVdt =-2k-(g-Ap V ■-/ p IT)2 . --- (1.26)

o o a

I n t e g r a t i o n of t h i s expression with r e s p e c t to time y i e l d s

R2-R 2 = 2k (g Ap V / p T 0 * t . (1.27)

\J O O 3

Thus the cloud area is a linear function of time. (1.27) can be written in dimensionless form as follows

where

(R/R l2 - 2k (t/t ) + 1 (1.28)

o o

t m R /U (1.29)

is a time scale and

Uo = (g Apo Ho / P ƒ * (1.30)

a velocity scale.

Thus R0 and U0 are "natural" scaling parameters for radial gravity spreading.

These scaling parameters can also be used to write the equation (1.24) for the front velocity in dimensionless form. The result for large times (t/tQ>>1) is:

U./U = (kt /2t)^ . (I.3D f o o

We see that the front velocity decreases inversely proportional to the square-root of time, which means that the cloud is decelerating.

The spreading "law" (1.31) cannot be valid for early times. Initially the cloud is at rest, so there must be a period in which the cloud accelerates. We will not deal with this acceleration phase in detail. It is useful however to give an estimate for its duration. This estimate is obtained as follows. Initially the static pressure force is used to accelerate both the cloud mass itself and an amount of "added air mass". We approximate the static pressure force which acts on the cloud edge by:

Fg = 2TTR H Pg (1.32)

and the inertial force by

- d Uf

Fi = P V d F - ( I'3 3 )

Assuming a balance between these two forces and using (1.21), we find that for time t=0

dU„ 2g Ap H 2p U

r f-i o ° a o ( ni | )

Ldt Jt=0 - p R ~ p t * U , i 4 ;

o o o o

This suggests that an appropriate time scale ta for the initial acceleration

phase is

t - ! 2 _ t . »-35)

Since normally p /2p i s of 0 ( 1 ) , t _ / t

_{o a a v}

i s of 0 ( 1 ) . A numerical example i l l u s

-t r a -t e s -the kind of -time s c a l e s we may expec-t for -the i n i -t a l a c c e l e r a -t i o n . Le-t

us take a cloud with i n i t i a l l y Ap /p = 1, H

= 10 m, R

= 10 m, then we find

U

= 10 m s , t

= 1 s and t

= 1 s. Thus the i n i t i a l a c c e l e r a t i o n phase i s

very short in t h i s example. A more d e t a i l e d discussion on t h i s subject i s given

in chapter I I I , section 6.3«

Next we proceed with the scaling of the mixing process. We have seen

before that observations show that the cloud edges are very d i s t i n c t . This

i n d i c a t e s t h a t l i t t l e mixing occurs through the cloud edge. On the other hand

the cloud top i s quite diffuse, which i n d i c a t e s t h a t s i g n i f i c a n t mixing may

occur t h e r e . Thus as a f i r s t appriximation we describe the mixing process by

dV/dt =

W , (1.36)

e

where We is an area-averaged top entrainment velocity to be estimated later.

Because V = irR2 H and dR/dt = Uf we may write (1.36) as

dH/dt = - 2(H/R) U + W . (1.37) f e

This shows t h a t the cloud height changes by two processes. The f i r s t term a t

the r i g h t hand s i d e gives the mean downward f l u i d motion caused by the

spreading and slumping process. The l a s t term shows t h a t the entrainment v e l ­

o c i t y i s equivalent to a mean upward motion of the cloud top r e l a t i v e to the

f l u i d motion. In section 1.1 we have seen t h a t the p o t e n t i a l energy of the

cloud i s proportional to the cloud height ( 1 . 1 2 ) . Thus the entrainment r e l a t i o n

(1.37 )" i s " d i r e c t l y 'linked with the energy "budget "of the cloud. ""To"'be' more p r e ­

c i s e , the energy budget of the cloud imposes c e r t a i n r e s t r i c t i o n s on the

scaling laws of W

. W

should not be modeled in such a way t h a t the p o t e n t i a l

energy increases above i t s i n i t i a l v a l u e . On the contrary, we expect that the

p o t e n t i a l energy and the cloud height decrease with time because of continuous

d i s s i p a t i o n of energy i n t o h e a t . Thus we expect that

W < 2(H/R) U_ . (1.38)

e f

Let us now look at the entrainment process i t s e l f . Entrainment occurs because

turbulence i s produced by the shear that i s r e l a t e d to the r a d i a l gravity

spreading. The strongest shears are present near the advancing leading edge.

Here turbulent eddies are created with an energy density of order \ p U * .

cl I

T h i s occurs a t a volume r a t e of o r d e r 2TTR H Uf. thus t h e e d d y - p r o d u c t i o n r a t e

should s c a l e as irp R H U 3 . This p r o d u c t i o n r a t e i s used p a r t l y for t h e

3. I

viscous dissipation of turbulent energy, partly for increasing the potential energy. From (1.12) and (1.37) it follows that the latter energy transformation occurs at a rate i g Ap V W . We now assume similarity in the sense that a

o o e

fixed fraction e of the turbulent energy production is used for increasing the potential energy. This leads to

X g A]T VQ Wg - e ir pa R H Uf3 . (1.39)

Using (1.11), (1.22) and (1.21) we easily obtain that

where

a ■ e k2 (1.41)

e

is an entrainment coefficient. Comparing (I.40) and (1.37) we see that the energy constraint is met for a < 1. Further it is interesting to see that in a mathematical sense the scaling law (I.40) for top entrainment is identical to the edge entrainment proposed e.g. by Van Ulden (1974). A major problem with (1.40) is that the value of a is not well established. Values for a ranging from 0.4 to 1.0 have been given in the literature. Minor variations in a have a significant effect on cloud concentrations for large times. This can be seen, when we substitute (1.40) for We in (1.36) and solve for V by using (1.24). The

result is

a

V/V = [2k(t/t ) + 1] e . (1.42)

o o

From this result we see that for great times V/VQ becomes very sensitive for

the value of a . Thus a significant amount of uncertainty exists here. More

discussion on this subject is given in chapter III. Despite the uncertainties in the scaling "laws" (1.28) and (1.42), these laws are useful for providing a first educated guess for the spreading and mixing of an isothermal dense cloud. We will illustrate this with an example.

Let us take instantaneous release of a cylindrical cloud with RQ = 10 m,

H„ = 10 m and Ap /p = 1 . Then U. = 10 m s" ' and t_ = 1 s. Let us assume that

o o a ° °

k = 1 and a = 0.5. With these values we find from (1.28) that after 50 s the e

cloud area has been increased by a factor 100, i.e. R is about 100 m. From (I.42) we see that the cloud volume has been increased by a factor 10. Thus average cloud concentration c and its relative density difference with air Ap/p are decreased to about 10? of their initial values. The cloud height is also decreased to about 10% of its initial value, i.e. to H = 1 m. By then, according to (1.24), the spreading velocity Uf = 1 m s"1. In atmospheric flow

conditions with low windspeeds this is a very significant velocity.

The evaluation after 20 minutes gives R = 500 m, H = 0.2 m, A~p7p = 0.02 and c = 0.02. Then the spreading velocity has decreased to LL. = 0.2 m s- 1,

which still is a significant spreading velocity.

From this numerical example we learn that a large dense cloud may spread over a very large area and becomes very shallow. At the same time significant self-mixing occurs. Also we see that the time scale of all this is quite small. From this it is apparent that gravity spreading and self-induced mixing are features that should be included in models for the dispersion of large dense clouds.

C h a p t e r I I

A DYNAMICAL INTEGRAL MODEL FOR TWO-DIMENSIONAL GRAVITY CURRENTS

11.1 Introduction

The gravity spreading phenomena due to negative buoyancy of dense clouds are not unique. Similar phenomena have been observed in other "man-made" flows and in a number of geophysical flows. Examples are:

- exchange flows in locks, which separate salt water and fresh water (e.g. Barr, 1967; Simpson and Britter, 1979);

- the spreading of oil on water (e.g. Abbott, 1961; Fay, 1969; Fanneldp and Waldman, 1972);

- atmospheric cold fronts (e.g. Schmidt, 1911).

For a recent review on gravity currents we refer to Simpson (1982).

Since gravity currents and especially two-dimensional gravity currents have been studied widely, it is useful to study the 2-D problem first before we enter into the details of the axisymmetric spreading problem. In this chapter we focus on horizontal bottom currents, which result from an instantaneous release of dense fluid. For reasons of simplicity we focus on shallow currents in which vertical accelerations can be neglected. We further assume that the ambient fluid is at rest and infinitely deep. Viscous forces are neglected.

The structure of this chapter is as follows. We start in section 2 with a review of studies made on "inertia-buoyancy currents" and discuss some general features of these currents. In section 3 we describe a new bulk model. This models consists of bulk equations for matter and momentum, which are derived from the basic equations of motion. In the sections 4 and 5 we present model results and an evaluation of these. In section 6 we give our conclusions.

11.2 History and problem definition

II.2.1 The leading edge of gravity currents

Major information on the conditions of the leading edge can be found in Schmidt (1911), Benjamin (1968), Simpson (1972) and Simpson and Britter (1979). In figure II.1 an arrested leading edge is shown.

t

t

figure II.1: The head of a steady gravity current (after Simpson and Britter, 1979).

Some characteristic features are:

1. At t h e leading edge a h e a d is present with a depth H-^ w h i c h is about twice the d e p t h Hn o f the c u r r e n t behind the h e a d .

2 . A n elevated forward s t a g n a t i o n p o i n t is p r e s e n t . B e l o w this point a n insig­ nificant flux Qi o f a m b i e n t f l u i d is e n t r a i n e d .

3 . Behind t h e h e a d a w a k e r e g i o n is p r e s e n t in which s i g n i f i c a n t mixing o c c u r s . 4 . In the h e a d a s i g n i f i c a n t internal f l o w is present. Near the surface d e n s e r

fluid m o v e s towards t h e h e a d with a velocity Uj|. In the w a k e region a m i x i n g layer m o v e s away from t h e h e a d .

S & B found that

U^.- .0....2..U., (II.2.1)

where U is the velocity of the current relative to the ambient fluid. 5. The dimensionless velocity Ch of the leading edge can be written as

Ch - U / / 8APhHn/Pa . (II.2.2)

where Ap. is the mean density difference behind the head, p the density of

h a the ambient fluid and Hn the densimetric mean depth behind the head. A

formal definition of the product of Ap and Hh is (Fay, 1980):

H,

APu H. Kh h

With t h i s d e f i n i t i o n the value of C

may be computed from experimental d a t a .

From Schmidt (1911), Benjamin's (1968) review, S & B (1979), Fannelefp e t a l .

(1980), Huppert and Simpson (1981) i t follows t h a t

Ch = 1 .15 ± 0.05, ( I I . 2 . 1 )

both for steady and unsteady gravity currents, provided the Reynolds number UH^/v > O(10^) and provided the current is deeply submerged.

II.2.2 Bulk properties of fixed volume releases

A major contribution to the understanding of the bulk properties of fixed volume releases has been given by Fanneldp and Waldman (1971, 1972) and later by Hoult (1972). Their aim was to describe the spreading of oil on water. They assumed (II.2.2) as a leading edge boundary condition and used the shallow-water equations to describe the interior of the current. They showed that similarity solutions to this set of equations exist for the layer averaged velocity u, for the local depth h and for the dimensionless velocity

C„ ■ U„//gAp"H/pQ , (II.2.5)

n I a

where Uf is the front velocity, Ap the mean density difference and H = V/X the

mean depth of the current. (V is the volume per unit width of current and X the length of the current.) Their results cannot be applied as such to the present problem, because entrainment and the internal flow were not included in their description. Also their solutions are only valid for large times, when the current has passed through its initial acceleration phase into the final deceleration phase. The problems of entrainment and initial acceleration have been dealt with - be it in a crude manner - by Van Ulden (1979). Van ULden assumed a rectangular shape with a linear velocity distribution and derived a bulk equation for dllf/dt. He included a static pressure force, a drag force and

an effective stress due to entrainment in this equation. However, objections may be raised against the way in which this equation was derived.

It is the purpose of this paper to improve on this. We will derive bulk continuity and momentum equations starting from the basic equations of motion. The now existing better understanding of the leading edge conditions will be employed in the model. Also some features of the approach by Fanneldp and

Waldman (1971 , 1972) will be included in the model. Furthermore a parameter­ ization for top entrainment will be presented.

II.3 A new bulk model for fixed volume releases

II.3.1 Introduction

In the remainder of this paper we deal with gravity currents which result from an instantaneous release of a volume V0 (per unit width) of a fluid with

density p that is greater than the density p of the ambient fluid by an

O cl

amount Ap . The two fluids are assumed incompressible and of equal temperature. The release is at the horizontal bottom at the beginning of an infinitely deep channel in which the ambient fluid is at rest. The initial volume has a length XQ and a mean depth H = V /X . After the release a gravity current develops of

the type that is shown in figure 2. Such currents develop when surface friction and side wall effects are negligible and the current is shallow. These condi­ tions are presupposed in this chapter. In the figure it is shown that a typical gravity current has a distinct head, followed by a long tail with a diffuse upper boundary. The local depth of the gravity current is defined by:

h(x,t) = 2 ƒ z Ap(x,z,t) dz / ƒ Ap(x,z,t) dz . (II.3.1)

Thus the local depth is twice the local densimetric mean depth z(x,t). For our model we will use an idealised representation of a real gravity current. Our idealised gravity current has a uniform density p and all dense material is below the local depth h(x,t). The volume per unit width of the current is given by

t

t

t

V = o/X h(x,t) dx (II.3.2)

and the horizontal momentum-integral per unit width is approximated by

M = o/X pc ü(x,t) h(x,t) dx , (II.3.3)

where ïï(x,t) is the layer averaged horizontal velocity. Equations for dV/dt and dM/dt form the basis of our dynamical integral model. Such equations will be derived in the sections II.3.2 and II.3-3. To obtain a full model closure also the integrals in (II.3.2) and (II.3.3) have to be evaluated. This is done in section II.3.4 for the tail region and in section II.3.5 for the head of the current. In these sections the volume and momentum of the tail and the head are expressed in terms of the length X of the current, the depth H^ in the origin, the depth Hn at the transition from tail to head and the velocity U = dX/dt of

the leading edge.

II.3.2 Equations for the volume of the current

As we have seen in II.2.1 little mixing occurs through the leading edge of the gravity current. On the other hand significant mixing occurs at the top of the current behind the head. Therefore we parameterize the mixing process as:

dV/dt = W X , (II.3.4) e

where We an effective entrainment velocity an X the area per unit width of the

current. Thus WeX can be interpreted as the volume flux of ambient fluid into

the gravity current.

Next we will derive a parameterisation for the entrainment velocity by considering the energy budget of the current. The potential energy of the current is given by

where

PE = 0/X Q/h ( x , t ) g Apc z dz dx (II.3.5)

Apc = pc - pa (II.3.6)

When h(x,t) does not vary too wildly with x the potential energy can be

approximated by

PE «..£ g Ap

V H , (II.3.7)

where

H = V/X (II.3.8)

is the average cloud depth.

As we have seen in Chapter I (eqn. (1.11)) the conservation of dense

material implies that

Ap

V = Ap

V

. (II.3.9)

Therefore

and

PE - i g Ap V H (II.3.10)

* o o

dPE/dt

-.\ g Ap V (dH/dt) . (II.3.11)

o o

Using (II.3.*») and (II.3.8) and using that

dX/dt = U

(II.3.12)

it is easily shown that

dH/dt =_^HU

/X + W

. _ (II.3.13)

Therefore

dPE/dt =

- \ g Ap V H U./X +.X g Ap^ V W . (II.3-1

*)

• ' o o f ' o o e

The physical interpretation of this result is that the potential energy

decreases due to slumping and increases due to entrainment. The latter process

corresponds with a transformation of turbulent energy into potential energy.

This transformation process is often called "the buoyant destruction of turbu­

lence". From (II.3.1H) we see that the entrainment velocity is known when the

buoyant destruction is known. To find the buoyant destruction rate we make the

closure assumption that buoyant destruction is proportional to the production of turbulent kinetic energy. This assumption has led to successful modeling of entrainment in the atmosphere and ocean (e.g. Tennekes and Driedonks, 1981). The production of turbulent energy in the present case can be derived from the analysis by Simpson and Britter (1979). They showed that shear production occurs mainly in the head region of the current, while destruction of eddies mainly occurs behind the head. Near the head, eddies are created with an energy density of o r d e r\ pa Uf . This occurs at a volume rate of order HnU. Thus the

production rate should scale as % p H, uj*. Using this we obtain for the

a h f buoyant destruction rate:

' X g Apo VQ We = \ e pa Hh Uf 3 , (II.3.15)

where e is a numerical coefficient. From this result we find with (II.3.'O and (II.3.9) for the entrainment rate:

where

dV/dt = E H , U./R1 (II.3.16) h f

Ri = gAp H/p U.2 (II.3.17)

C cl X

is a bulk Richardson number and where e is an empirical coefficient, this re­ sult resembles the traditional scaling of side entrainment (Van Ulden, 197^; Fay, 1980) which reads or the two-dimensional case

dV/dt = a HU„ , (II.3.18) e f

where a is a constant. The physical meaning of our result is, however, completely different. In our model only the production of turbulent kinetic energy occurs at the leading edge, but the subsequent entrainment occurs at the top of the current. Thus turbulent eddies are mainly created near the leading edge. While travelling away from the leading edge they lose their kinetic energy and increase the potential energy of the current.

The estimation of our entrainment coefficient e is not easy. No data seem to be available for thé 2-dimensional case. As a preliminary value we will use e = 0.6. This value lies within the range of entrainment coefficients found for the axisymmetric case (see chapter III for a full discussion on this subject).

II.3.3 The momentum-integral equation

In this section we will derive an equation for the horizontal

momentum-integral of the gravity current by considering the momentum budget in a time

varying control volume, which just encloses the gravity current. This budget

reads:

dM/dt = -

/

/

(3P/3x) dx dz . (II.3.19)

In this equation P is the pressure, while X+ and h+ denote, that the integra­

tions are extended to include the outer boundary of the current. In this

momentum budget it is presupposed that shear-stresses and momentum fluxes

vanish at the boundaries of the control volume. In this section we will

evaluate the right hand side of this momentum budget, which represents the

pressure forces that act on the current.

It is convenient for this evaluation to introduce the pressure disturb­

ance II, which is defined as:

n(x,z,t)

= P(x,z,t) - P j z ) , (II.3.20)

where P (z) is the pressure of the undisturbed ambient fluid at the height z.

Since 3P /3x = 0, it follows that

- ;

;

(3P/3X)

dx dz = - ;

;

on/ax) dx dz . ( n . 3 . 2 1 )

0 0 0 0

Since 3h+/3x and II at the upper boundary are small in the bulk of the current

we may write:

X+h+ X+ 3 h+ h+(o,t) h+(X+,t)

-I I (3Il/3x)dxdz = -ƒ — ( ƒ ndz)dx = ƒ n(0,z,t) -ƒ n(X+,z,t)dz

00 o

o o o (n.3.22)

Thus the pressure force on the current is determined by the pressure disturb­

ances in the origin and at the leading edge.

Three different processes produce these pressure disturbances:

- the negative buoyancy of the dense fluid gives rise to an increase II in the

static pressure in the current and a corresponding static pressure force F

;

- the advancing head of the current intrudes into stagnant ambient fluid. This

produces a dynamical pressure n. in front of the head and a corresponding

dynamical drag force Fd;

- the acceleration of the current applies momentum changes in the ambient fluid. This leads to an acceleration reaction pressure n in front of the head and a corresponding reaction force Fa.

In the following we will describe these three pressure forces.

The static pressure disturbance vanishes in X+, because no dense fluid is present at this distance from the origin. Thus the static pressure force is given by:

F

s

/

n

(o,z,t) dz . (ii.3.23)

Since n obeys the hydrostatic relation

S *

311 /3z = - gAp (II.3.24) s c

it can be shown that the static pressure force equals:

Fs =\ gApc Hfc2 , (II.3.25)

where Hfc = h(o,t) is the current depth in the origin (see fig. II.2).

The dynamic pressure vanishes in the origin. Thus the dynamic drag force is determined by the flow conditions near the leading edge. The dynamic

2

pressure difference over the leading edge scales as \ p U_ , while the scale

cL I

depth for the leading edge is Hh (see fig. II.2). Therefore we parameterize the

dynamic drag force as:

Fd = "-* cf pa Uf2 Hh ' ' (II.3.26)

where. cf is a drag coefficient of 0(1) which will be estimated in section

II.3.5.

The acceleration reaction force is written as:

Fa = - dMy/dt , (II.3.27)

where My is the virtual momentum of the ambient fluid (Batchelor, 1967, 6.4).

To parameterize this virtual momentum we use the resemblance between our gravity current and an elliptical cylinder with aspect ratio H/X. For such a

2

cylinder the virtual momentum equals p H U (Batchelor, 1967, 6.6). Therefore

3. X

we parameterize My as:

My = a pa H2 Uf , (II.3.28).

where a is an empirical coefficient of 0(1). The exact value of this

coefficient is not important for a shallow gravity current, because Mv is small

in comparison with the total momentum M of the current. In this study we use a = 2.

We are now ready for our momentum-integral equation. Using that

dM/dt = Fs + Fd + Fa we easily find from (II.3.25), (II.3.26) and (II.3.27)

that:

d(M+My)/dt *>.\ gApcHt2 -X cfPaHhUf 2 . (II.3.29)

This equation and the continuity equation (II.3.16) are the bulk rate equations for the total gravity current. In these equations the volume integral (II.3.2) and the momentum integral (II.3-3) still have to be specified. We will do so in

the following sections. ^

II.3.^ The horizontal distributions of layer depth and layer averaged velocity in the tail of the current

It is the purpose of this section to evaluate for the tail of the current the volume-integral

_._.h

Vt(t) ■ •ƒ h(x,t) dx (II.3.30)

o

and the momentum-integral

Xh

M ( t ) B .ƒ p u ( x , t ) h ( x , t ) dx ( I I . 3 . 3 D

o

In particular we want to express these integrals in terms of the model

variables H^, Hh, Xn and Un. In order to do so we need approximations to the

functions h(x,t) and u(x,t). During the early development of the current these are difficult to obtain, but quite soon the current becomes shallow enough that

the shallow layer equations are applicable to the flow in the tail of the current. In the present problem the following equations apply:

and

DÜ

g

3u

3x

c

c

w

h

3h

c

( I I . 3 . 3 2 )

X 6 ■ % • ( I I . 3 . 3 3 )

In these equations D/Dt = 3/3t + u 3/3x, w

i s the l o c a l entrainment velocity

and X 6 3u /3x a momentum flux gradient term that accounts for the fact t h a t

the horizontal velocity varies with height. 6 i s an empirical constant to be

estimated l a t e r . I t should be noted t h a t ( I I . 3 . 3 3 ) i s f u l l y c o n s i s t e n t with our

bulk equation ( I I . 3 . 2 9 ) , provided the current i s shallow. Also in ( I I . 3 . 2 9 ) an

e f f e c t i v e s t r e s s gradient term s i m i l a r to t h a t in ( I I . 3 . 3 3 ) is hidden in the

dM/dt term. This can be checked by evaluating the time d e r i v a t i v e of ( I I . 3 . 2 0 ) .

The physical meaning of the d e c e l e r a t i n g s t r e s s gradient term i s simply t h i s .

Entrainment does not affect the t o t a l momentum - since no momentum i s

entrained - , but i t increase the t o t a l mass of the c u r r e n t . This n e c e s s a r i l y

causes a decrease in the mean v e l o c i t y . Thus entrainment leads t o a d e c e l e r a ­

t i o n term in the equation for the mean v e l o c i t y . We will now derive approximate

solutions to ( I I . 3 . 3 2 ) and ( I I . 3 . 3 3 ) by making two s i m i l a r i t y assumptions. The

f i r s t i s that the shape of the c u r r e n t i s quasi-conserved in time i . e . t h a t

I B . l d H ( I I 3 3H)

h Dt H dt ' U l . J . J i ;

virtually independent from x. The second assumption is that the layer averaged density difference remains horizontally uniform and equal to Ap . This requires that

w /h = W /H , (II.3.35) e e

v i r t u a l l y independent from x (see also section I I . 3 . 1 ) . I t then follows t h a t

3u/3x = U./X. ( I I . 3 . 3 6 )

h h

and that

ü = x U./X. ( I I . 3 . 3 7 )

h h

I t also follows from ( I I . 3 . 3 3 ) - ( H . 3 . 3 7 ) t h a t 3h/3x i s a l i n e a r function of x

t h a t vanishes in x = 0. Using the boundary conditions h = H

for x = 0 and

h = Ht + (Hh-Hfc) (x/Xh)2 (II.3.38)

The solutions for ü and h happen to be of the same form as those obtained by Fanneldp and Waldman (1971, 1972). However there are two differences. In our model Hj. and Hn are independent variables that are determined by the dynamics

of the gravity current. Further we use (II.3.37) and (II.3.38) only to estimate the volume and momentum integrals of the tail. From (II.3.30), (II.3.31) and (II.3.37), (II.3.38) we easily obtain

and

Vfc = j (2 Ht + Hh) Xh (II.3.39)

M

t " ï

c

t

V

h

h

-

^

This completes our description of the tail.

II.3.5 The head of the current

The shallow layer equations are not applicable to the head of the

current. Instead we use the momentum-integral approach that we applied to the total current (II.3.19). The force balance for the head looks as follows. The static pressure force follows from the integration of the static pressure gradient over the head region and equals:

2 s '2 °"Kc "h

Fo -X gAp„ H ^ (II.3.41 )

The dynamic pressure force is

Fd " --* cf Pa Hh U/ * ( I I'3'1 4 2 )

Furthermore there is a momentum flux into the head due to the internal current in the head (figure II.1). Near the surface the inward flow Ujj carries positive momentum into the head. The return flow IU carries negative momentum out the head. So the net effect of the internal flow is a positive momentum flux into

the head. Assuming U = U. and h. = \ H we find that this flux is crudely

Using ü. =» 0.2 Uf (II.2.1) we may write this as

Qh " % ö pa Hh Uf2 (II.3.44)

where

6 - 0.08 . (II.3.45)

Thus 6 is an empirical coefficient which characterizes the non-uniformity of the vertical velocity profile. It has the same meaning as in (II.3.37).

We neglect the inertial terms in the momentum-integral equation for the head. It can be shown that these terms are normally small in comparison with

the other terms. Thus we assume that the head is in a quasi-steady state. This assumption is supported by experiment. We have seen in section II.2.1 that the dimensionless leading edge velocity Cn = 1.15 ± 0.05 both for steady and for

unsteady currents (II.2.4). Our momentum-integral equation now reads:

0 - \ gApcHh2 -X cfPaHhUf 2 -X öPaHhUf 2 (II.3. 46)

It follows from this equation that the dimensionless velocity Cn defined by

(II.2.2) is given by

C - 1 / / C.-Ó . (II.3.47) n f

Since Cn and 6 are known, the value of c^ can be estimated from this equation.

The result-is

cf = 0.84 + 0.07 (II.3.48)

With (II.3.44)-(II.3.48) we have specified the important dynamical leading-edge boundary conditions.

We conclude this section with the specification of the volume and momen­ tum of the head. We allow the volume of the head to vary in time, but assume that its shape remains unchanged. In section II.2 we have seen that the depth H-^ of the head is about twice the depth Hn behind it. We also assume that its

length scales with Hn. Thus we write

Experimental data suggest that a = 2. The volume of the head now is written as

Vh - b Hh 2 , (II.3.50)

where b = 4. To estimate the momentum of the head we assume that (as in the tail) the layer averaged velocity increases linearly with x. It then follows that

Un = UfXn/X . (II.3.51)

Using ( I I . 3 . 4 9 ) - ( I I . 3 . 5 1 ) we f i n d for the momentum of t h e head. 2

'c "f"h v" 'z " " h '

M^ = bp„ U ^ (X - \ a H J / X . ( I I . 3 . 5 2 )

This completes our description of the head. Together, the bulk equations derived in II.3.2 and II.3.3 and the equations for the tail and the head

derived in II.3.4 and II.3.5 form a closed set. We will summarize the final set of equations in the next section.

II.3.6 Final model equations

In the final model equations we use the dimensionless density difference A = Ap /p and the velocity integral M = (Mt+Mh+My)/p . The model has 7

variables, i.e. X, V, M , Uf, Hn, Hfc and A and 7 equations, namely'3 rate

equations and H diagnostic equations. These follow now. The first rate equation

results from the definition of Uj. and reads: _ _

dX/dt - Uf . (II.3.53)

The second rate equation follows from (11.3.^) and (II.3.15):

dV/dt - e H.XU_3/gA* V , (II.3.51»)

n f o o

where e = 0 . 6 i s an e n t r a i n m e n t c o e f f i c i e n t , and A and V a r e t h e i n i t i a l

* ° °

v a l u e s of A and V. The t h i r d r a t e e q u a t i o n f o l l o w s from ( I I . 3 . 2 9 ) :