-------
FOREWORD
Since 1974 the EPA Office of Radiation Programs (ORP) has
initiated studies, pertaining to ocean disposal of low-level
radioactive wastes (LLW), in response to the Marine Protection,
Research and Sanctuaries Act (Public Law 92-532) of 1972. The
studies were conducted to: determine the condition of LLW
containers in the ocean (after disposal on the seafloor and
prolonged exposure to the marine environment); assess any effects
to man or the environment from those disposals; and obtain the
information/data needed to promulgate regulations and criteria
for any future disposals.
This report discusses and summarizes the results of an ORP
study completed under an interagency agreement with the
Department of Energy (DOE). The DOE Battelle Pacific Northwest
Laboratory (BPNL) was tasked to: (1) identify regional ocean
models that could be used to determine effects from ocean
disposals of LLW; (2) evaluate the mathematical representations
for eddy viscosity/diffusivity, and other, coefficients in those
models; (3) determine the applicability of using a k-£
turbulence model; and, (4) assess the feasibility of applying one
of the identified models to the mid-Atlantic Ocean, 2800 m and
3800 m sites that were used previously for disposal of LLW.
The Agency invites all readers of this report to send
comments or suggestions to Mr. Martin P. Halper, Director,
Analysis and Support Division (ANR-461), Office of Radiation
Programs, Environmental Protection Agency, Washington, DC 20460.
Richard^J. Guimond, Director
Office of Radiation Programs
111
-------
SUMMARY
Pacific Northwest Laboratory (PNL) performed a study for the
U.S. Environmental Protection Agency's (EPA) Office of Radiation
Programs (ORP) to: (1) identify candidate models for regional
modeling of low-level radioactive waste (LLW) ocean disposal sites
in the mid-Atlantic ocean; (2) evaluate mathematical representation
of the models' eddy viscosity/dispersion coefficients; and, (3)
evaluate the adequacy of the k-fc turbulence model and the
feasibility of one of the candidate models, TEMPEST•/FLESCOT*f to
deep ocean applications on a preliminary basis.
PNL identified the TEMPEST/FLESCOT, FLOWER, Blumberg's, and
RMA 10 models as appropriate candidates for the regional
radionuclide modeling. Among these modesls, TEMPEST/FLESCOT is
currently the only model that solves distributions of flow,
turbulence (with the k-e model), salinity, water temperature,
sediment, disssolved contaminants, and sediment-sorbed
contaminants.
Solving the Navier-Stokes equations using higher order
correlations is not practical for regional modeling because of the
prohibitive computational requirements; therefore, the turbulence
modeling approach is more practical.
PNL applied the three-dimensional code, TEMPEST/FLESCOT with
the k-€ model, to a very simple, hypothetical, two-dimensional,
deep-ocean case, producing at least qualitatively appropriate
results. However, more detailed testing should be performed to
further test the code.
Copyright by Battelle Memorial Institute, Columbus, Ohio.
-------
CONTENTS
FOREWORD ill
SUMMARY V
LIST OF FIGURES viii
LIST OF TABLES xi
1. INTRODUCTION 1.1
2. REGIONAL MODELS FOR LOW-LEVEL RADIOACTIVE WASTE
SITES IN THE MID-ATLANTIC 2.1
3. REVIEW AND EVALUATION OF MATHEMATICAL REPRESENTATION
FOR VISCOSITY AND DISPERSION COEFFICIENTS 3.1
3.1 ZERO-EQUATION MODELS 3.3
3.2 ONE-EQUATION MODELS 3.8
3.3 TWO-EQUATION MODELS 3.11
4. MODEL DESCRIPTION 4.1
4.1 TEMPEST/FLESCOT MODEL 4.1
4.2 HYDRODYNAMIC EQUATIONS 4.3
4.3 HEAT TRANSFER 4.7
4.4 MASS TRANSPORT 4.8
4.5 TURBULENCE EQUATIONS 4.8
4.6 SEDIMENT TRANSPORT EQUATION 4.10
4.7 DISSOLVED CONTAMINANT TRANSPORT EQUATION 4.12
4.8 SEDIMENT-SORBED CONTAMINANT TRANSPORT EQUATION 4.13
5. MODEL APPLICATIONS 5.1
5.1 EDDY VISCOSITY DISTRIBUTIONS 5.1
5.1.1 Nonstratified Case (Case 1) 5.1
5.1.2 Stratified Cases (Cases 2 and 3) 5.3
5.2 RADIONUCLIDE"DISTRIBUTIONS 5.6
5.2.1 Simulation Conditions 5.6
5.2.2 Simulation Results of Stratified Case 4 5.7
5.2.3 Simulation Results of Stratified Case 5 5.9
5.2.4 Simulation Results of Stratified Case 6 5.10
6. SUMMARY AND CONCLUSIONS 6.1
7. REFERENCES 7.1
vii
-------
FIGURES
No. Page
5-1 Calculated Turbulent Kinetic Energy Over the
Whole Flow Region for Nonstratified Case l 5.11
5-2 Calculated Turbulent Energy Dissipation Over
the Whole Flow Region for Nonstratified Case 1 5.12
5-3 Calculated Eddy Viscosity Over the Whole Flow
Region for Nonstratified Case 1 5.13
5-4 Calculated Turbulent Kinetic Energy Near
the Bottom for Nonstratified Case 1 5.14
5-5 Calculated Turbulent Energy Dissipation Near
the Bottom for Nonstratified Case 1 5.15
5-6 Calculated Eddy Viscosity Near the Bottom
for Nonstratified Case 1 5.16
5-7 Vertical Distribution of Calculated Turbulent
Kinetic Energy Near the Bottom for
Nonstratified Case 1 5.17
5-8 Vertical Distribution of Calculated Turbulent
Energy Dissipation Near the Bottom for
Nonstratified Case 1 5.18
5-9 Vertical Eddy Viscosity Near the Bottom
for Nonstratified Case 1 5.19
5-10 Calculated Velocity Distribution Near the
Bottom for Nonstratified Case 1 5.20
5-11 Calculated Turbulent Kinetic Energy Over
the Whole Flow Region for Stratified Case 2 5.21
5-12 Calculated Turbulent Energy Dissipation Over
the Whole Flow Region for Stratified Case 2 5.22
5-13 Calculated Eddy Viscosity Over the Whole
Flow Region for Stratified Case 2 5.23
5-14 Calculated Turbulent Kinetic Energy Near the
Bottom for Stratified Case 2 5.24
viii
-------
FIGURES (Continued)
No. Page
5-15 Calculated Turbulent Energy Dissipation
Near the Bottom for Stratified Case 2 5.25
5-16 Calculated Eddy Viscosity Near the Bottom
for Stratified Case 2 5.26
5-17 Vertical Distribution of Calculated Turbulent
Kinetic Energy Near the Bottom for
Stratified Case 2 5.27
5-18 Vertical Distribution of Calculated Turbulent
Energy Dissipation Near the Bottom for
Stratified Case 2 5.28
5-19 Vertical Eddy Viscosity Near the Bottom for
Stratified Case 2 5.29
5-20 Calculated Velocity Distribution Near the
Bottom for Stratified Case 2 5.30
5-21 Calculated Turbulent Kinetic Energy Over the
Whole Flow Region for Stratified Case 3 5.31
5-22 Calculated Turbulent Energy Dissipation Over
the Whole Flow Region for Stratified Case 3 5.32
5-23 Calculated Eddy Viscosity Over the Whole Flow
Region for Stratified Case 3 5.33
5-24 Calculated ^Turbulent Kinetic Energy Near the
Bottom for Stratified Case 3 5.34
5-25 Calculated Turbulent Energy Dissipation Near
the Bottom for Stratified Case 3 5.35
5-26 Calculated Eddy Viscosity Near the Bottom
for Stratified Case 3 5.36
5-27 Vertical Distribution of Calculated Turbulent
Kinetic Energy Near the Bottom for
Stratified Case 3 5.37
IX
-------
FIGURES (Continued)
No. Page
5-28 Vertical Distribution of Calculated Turbulent
Energy Dissipation Near the Bottom for
Stratified Case 3 5.38
5-29 Vertical Eddy Viscosity Near the Bottom for
Stratified Case 3 5.39
•lji
5-30 Calculated Velocity Distribution Near the
Bottom for Stratified Case 3 5.40
5-31 Calculated Dissolved Radionuclide Concentration
(in pCi/L) Near the Bottom for Stratified Case 4 5.41
5-32 Calculated Suspended-Silt-Sorbed Radionuclide
Concentration (in pCi/g) Near the Bottom for
Stratified Case 4 5.42
5-33 Calculated Suspended-Clay-Sorbed Radionuclide
Concentration (in pCi/g) Near the Bottom for
Stratified Case 4 5.43
5-34 Calculated Dissolved Radionuclide Concentration
(in pCi/L) Near the Bottom for Stratified Case 5 5.44
5-35 Calculated Suspended-Silt-Sorbed Radionuclide
Concentration (in pCi/g) Near the Bottom for
Stratified Case 5 5.45
5-36 Calculated Suspended-Clay-Sorbed Radionuclide
Concentration (in pCi/g) Near the Bottom for
Stratified Case 5' 5.46
5-37 Calculated Dissolved Radionuclide Concentration
(in pCi/L) Near the Bottom for Stratified Case 6 5.47
5-38 Calculated Suspended-Silt-Sorbed Radionuclide
"concentration (in pCi/g) Near the Bottom for
Stratified Case 6 5.48
5-39 Calculated Suspended-Clay-Sorbed Radionuclide
Concentration (in pCi/g) Near the Bottom for
Stratified Case 6 5.49
x
-------
TABLES
No. Page
2-1 Summary of Reviewed Hydrodynamic
and Transport Models 2-4
3-1 Variation in Eddy Diffusivity 3-4
5-1 Initial Water Temperature and Salinity
Distribution for Case 2 5-3
5-2 Initial Water Temperature and Salinity
Distribution for Case 3 5-4
5-3 Simulation Conditions and Model Parameters
for Cases 4 through 6 5-8
5-4 Calculated Radionuclide Concentration Sorbed
by Sediment in the Top 3-cm Bed Layer
after 6 h Simulation 5-9
-------
1. INTRODUCTION
The U.S. Environmental Protection Agency (EPA) is investigate
ing the potential environmental impacts and health risks from low-
level radioactive waste (LLW) disposal on the deep-ocean bottom in
the mid-Atlantic. Mathematical models can assist in predicting
migration and fate of radionuclides from the LLW disposal site.
Pacific Northwest Laboratory (PNL) performed this study for
EPA with the following objectives: (1) identify candidate models
useful for regional modeling of the LLW ocean disposal sites; (2)
review and evaluate mathematical representations of eddy viscosity
and dispersion coefficients; and, (3) apply the TEMPEST/FLESCOT
model to a very simple hypothetical deep-ocean disposal case to
evaluate the adequacy of the k-e turbulent model on a preliminary
basis and the feasibility of the TEMPEST/FLESCOT model to the deep-
ocean applications.
Section 2 of this report summarizes PNL's review of represen-
tative mathematical models, Section 3 reviews and evaluates the
mathematical equations for viscosity and dispersion coefficients,
and Section 4 gives a description of the TEMPEST/FLESCOT model and
its governing equations. Section 5 describes the applications of
the TEMPEST/FLESCOT model in six cases. Conclusions are given in
Section 6.
1.1
-------
2. REGIONAL MODELS FOR LOW-LEVEL RADIOACTIVE WASTE
SITES IN THE MID-ATLANTIC
Low-level radioactive wastes (LLW) released from the
2800-m or 3800-m sites in the mid-Atlantic Ocean will not be
distributed uniformly in either the vertical or horizontal
direction on local (e.g., within several kilometers square) and
regional scales (e.g., up to several hundred kilometers square).
The physical-oceanographic features also exhibit complex
three-dimensional behavior. Therefore, as discussed in Onishi et
al. (1987), three-dimensional coupled hydrodynamic-mass/energy
transport models with proper representations of oceanic eddy
viscosity and dispersion processes are required as regional
radionuclide transport models.
The Pacific Northwest Laboratory reviewed some
representative mathematical models to evaluate their
applicabilities and limitations for the LLW ocean disposal assess-
ment (Onishi et al. 1987). Table 2-1 shows the review summary of
these models. As shown in Table 2-1, with suitable model
modifications or inclusions, such as turbulence closure, enhanced
sediment transport, radionuclide, and/or curvilinear coordinate
system setup, the FLESCOT model (Onishi and Trent 1982), the FLOWER
model (Eraslan et al. 1983), and Blumberg's model (Blumberg and
Herring 1986) would be appropriate candidates for regional radio-
nuclide modeling to predict the transport and dispersion of LLW
disposed in the 2800-m and 3800-m ocean sites (Onishi et al. 1987).
Although the RMA 10 model (King 1982) does not incorporate a
turbulence closure scheme, this model, with some modifications, is
also an appropriate candidate for regional radionuclide modeling.
These models are also applicable to a local scale covering a few
square kilometers, such as the area within the 2800-m and 3800-m
sites. For addressing the long-term, global-scale impact of LLW
2.1
-------
disposed in these sites, the output of these regional models is to
be provided to global-scale models as boundary conditions or
sink/source terms appear to be the only ones which have
comprehensive descriptions of transport and fate of sediments,
dissolved contaminants, and sediment-sorbed contaminants.
Among the regional-scale models shown in Table 2-1,
FLESCOT solves distributions of flow, turbulence, salinity, water
temperature, sediments, dissolved contaminants (e.g.,
radionuclides, hazardous chemicals), and sediment-sorbed
contaminants. The sediment and contaminant transport is simulated
in both ocean water column and bottom sediments. Thus, the
FLESCOT model can be a candidate model to be applied to the 2800
and 3800-m sites to predict the transport and accumulation of LLW
on regional and local scales and to address the criteria (B) 1,
(D) 2, and (J) as 3 listed in Section 424(i)(l) of Public Law
97-424.
The main strengths of the FLESCOT model for the deep-
ocean LLW disposal site assessment are listed here. The code:
• is time-dependent, three-dimensional
• calculates hydrodynamie and energy transport
• calculates sediment transport for both cohesive and noncohesive
sediment
'' An analysis of the environmental impact of the disposed
action, at the site at which the applicant desires to dispose
of material, upon human health and welfare and marine life.
\
*® An analysis of the resulting environmental and economic
conditions if the containers fail to contain the radioactive
waste material when initially disposed at the specific site.
^ A comprehensive monitoring plan to be carried out by the
applicant to determine the full effect of the disposal on the
marine environment, living resources or human health, etc.
2.2
-------
• has the k-e turbulence model, as a turbulent closure scheme,
which needs to be tested in this study
simulates both dissolved and sediment-sorbed radionuclide
transport with interactions between radionuclides and
sediments
• predicts sediment and radionuclide distributions at the ocean
bottom with ocean bottom data
• has been applied to the coastal and estuarine environment to
calculate flow and mass/energy transport, including sediment
and contaminant transport
is computationally efficient through the extensive use of
vectorizing.
The major limitation of the FLESCOT code is the lack of model
applications to the deep-ocean environment. The code must be
tested to examine suitability of representations of turbulence/eddy
viscosity and dispersion in the deep-ocean (such as the 2800- and
3800-m sites) environment.
-------
TABLE 2.1. Summary of Reviewed Hydrodynamic and Transport Models (Onishi et al. 1987)
Simulated Substances Sediment-
Model
Simons (1973)
Leenderste and
Liu (1975)
Ueatherly and
Martin (1978)
RA10
(King 1982)
FLESCOT©
(Onishi and
Trent 1982)
TEMPEST©
(Trent
et al. 1989)
FLOWER
(Eras I an
et al. 1983)
REHIXCS
(Freitas
et al. 1985).
Btunberg and
Herring (1986)
Sheng et al.
(1978)
Da vies (1980)
Tee (1981)
Turbulent Water Dissolved Particulate Contaminant Time
Flow Energy Salinity Temperature Sediment Contaminant Contaminant Interactions Dependency(a)
XX U
XXX U
XX X U 3 R
X X X X X X U
XXXXXX X X U
U
X X X X U
X X X X U
XX U
X X X X U
X U
X U
X U
Solution
Dimensions(b)
3
3
FD
3
3
3
3
3
3
3
3
3
3
Scale(c)
R
R
R
R
R
R
R
R
R
R
R
R
Technique(d)
FD
FD
FE
FD
FD
FD
DE
FD
FD
FD
FE
FD/HC
(a) "S" and "U" denote steady and unsteady states, respectively.
(b) "2", "3", and "C" denote two-, three-dimensional, and compartment models, respectively.
(c) "R" and "G" denote regional/local and global scales, respectively.
(d) "FD", "FE", "MC", "I", "A", and "OE" denote finite differences, finite element, method of characteristics, integration method, analytical,
and discrete element models, respectively.
-------
TABLE 2.1. (contd)
Simulated Substances
Sediment-
Model
CAFE
(Wang and
Conner 1975)
RMA2
(Norton and
King 1977)
FETRA
(Onishi 1981)
Spaulding and
Ravish (1984)
Walker
et at. (1985)
EXAMS
(Smith et al .
1977)
WASP
(Ditoro et al.
1981)
Bryan (1969)
Holland and
Liu (1975)
Webb and
Morley (1973)
Shepard (1978)
MARINRAD
(Koptik
et al. 1984)
MARKA
(Robinson and
Marietta 1985)
Turbulent Water Dissolved Particulate Contaminant Time
__ Flow Energy Salinity Temperature Sediment Contaminant Contaminant Interactions Dependency(a)
X U
X U
X X X X U
XX U
X U
XXX S
X X X X U
XXX U
X U
X U
X S
X X X U
X X X U
Solution
Dimensions(b)
2
2
2
3
2
C
c
3
2
3
3
C
C
Scale(c)
R
R
R
R
R
R
R
G
G
G
G
G
G
Technique(d)
FE
FE
FE
AFD
A
I
I
FD
FD
A
A
I
I
-------
3. REVIEW AND EVALUATION OF MATHEMATICAL REPRESENTATION FOR
VISCOSITY AND DISPERSION COEFFICIENTS
Determination of flow and constituent concentrations requires
the solution of the governing equations for continuity, momentum,
and scalar transport (temperature, salinity, sediments, and
radionuclides) shown below. To solve these equations, the terms
u.u. and u.0, known as the Reynolds stresses or turbulence
correlations, must be determined.
<3U.
Continuity ^T = ° (3.1)
i i i M
Momentum f^ + Uj -L + f. = JL |L_ (3>2)
J • 1
J J r
Scalar Transport
(temperature, concentration) 4£ + U. |^- = -£- (^- - uT0) + S (33)
dt i ax. ax1 \ dx. i 0
where U = mean flow velocity
F. = Coriolis parameter
p = density
pr - reference density
P = instantaneous static pressure
v = molecular viscosity
u = turbulent fluctuating velocity
g, = gravity in the x, direction
0 = scalar quantity
A = molecular diffusivity
S, = volumetric source term.
0
3.1
-------
Two methods are available to determine u( uj and U| 0. One
method for determining turbulence correlations uses exact
equations. In this method, the exact equations create higher order
correlations that must then be solved or modeled. Mellor (1973)
developed a method for modeling the second-order turbulence
correlations using the exact equations and modeling higher order
terms. Mellor and Yamada (1974, 1982) developed a series of models
using the same method. However, these higher order methods require
extensive computational time and are not very practical for
modeling with a large number of computational cells.
The second method, which probably is simpler, is to use the
concept of eddy viscosity/diffusivity to define the turbulence
correlations. The concept of eddy viscosity was developed by
assuming that the eddy viscosity is similar to molecular viscosity,
but eddy viscosity is based on the state of the turbulence and not
the properties of the fluid. Rodi (1984) developed the concept of
eddy viscosity and diffusivity more fully. The simplest method is
to use values specified directly for the viscosity and assume that
the diffusivity is a function of the eddy viscosity and the
turbulent Prandtl number for heat transport or the turbulent
Schmidt number for mass transport (Rodi 1984)
(3-4)
where k = eddy diffusivity
E = eddy viscosity
a = turbulent Prandtl or Schmidt number.
However, it is not reasonable to have a constant eddy
viscosity in a deep ocean. Thus, we will discuss how to obtain
nonhomogeneous , anisotropic eddy viscosity/diffusivity. Three
3.2
-------
types of models use the eddy viscosity/diffusivity concept to
attempt to predict the values of the turbulence correlations. One
way to specify the eddy viscosity is to use a functional
relationship based on depths and shear velocities similar to
that developed by Munk and Anderson (1948). Because this method
does not use a partial difference equation to express turbulence
process, it is called a zero-equation model. This relationship can
be modified to take into account the stabilizing effects of
stratification.
One-equation models are the second way of determining eddy
viscosity. One-equation models work well where the length scale
can be clearly defined, such as in the benthic boundary layer. If
the length scale is not as clearly defined, problems arise, and
thus the third meth'od, two-equation models, was developed. The k-
e model is a two-equation model developed to solve production,
dissipation, advection, and diffusion of turbulent kinetic energy,
which are then related to the turbulence correlations.
3.1 ZERO-EQUATION MODELS
Equations describing the mean velocity field can be developed
from the Navier-Stokes equation. Solving these equations requires
that values be obtained for the turbulence correlations, uj~~uj and
Uj
where k = turbulence kinetic energy.
3.3
-------
All terms are known except the eddy viscosity- The first method
proposed here is zero-equation models, which means they do not
require the solution of any additional partial differential
equations to resolve the flow field. The use of zero-equation
models is limited by the fact that the transport of turbulence is
not taken into account.
Specifying the eddy viscosity is one method used to obtain
a closure for the turbulence equations. This can be done by
specifying a constant eddy viscosity throughout the flow field or
specifying the values for specific regions in the flow field.
Specifying a constant-eddy viscosity is common, but it is not a
true turbulence model, and the results are quite crude. In
applications in the deep ocean, this method does not give an
accurate representation of what is actually occurring.
As seen in Table 3-1, the eddy diffusivity varies by several
orders of magnitude over the depth of the ocean, and because eddy
diffusivity is directly related to the eddy viscosity, one can
expect the eddy viscosity to vary also. The Prandtl or Schmidt
number is assumed to be constant in most cases. The specification
of regional values is an improvement, but the result is still
rather crude.
TABLE 3-1. Variation in Eddy Diffusivity
(U.S. Atomic Energy Commission 1975)
Layer
Surface
Intermediate
Deep
Depth, m
75
500
2000
Vertical Eddy
Diffusivity, Kv, cm2/s
50.0
0.1
1.0
3.4
-------
An alternative to specifying values of eddy
viscosity/diffusivity directly is the development of a functional
relationship for the variance of the eddy viscosity/diffusivity
throughout the depth. This concept has been used in estuaries and
shallow ocean regions, but no study has been done for the deep
oceans. In shallow regions where studies have been done, the flow
never develops beyond the boundary layer, whether it is on the
surface or the bottom. This results in the conclusion that it
should be possible to use the current functional relationships for
the boundary layers; however, these relationships will not apply
in the interior. A new relationship should be developed;
alternatively, it might be possible to use a constant value of eddy
viscosity/diffusivity, because not much variation is expected in
the interior regions.
c
Many methods have been developed to determine the
value for the vertical eddy diffusivity for an unstratified
boundary layer. Because stratification can play a major role in
the stabilization of turbulent flows, methods have also been
developed for modifying the eddy diffusivity as a function of the
Richardson number (Ri). One of the first such methods was
developed by Munk and Anderson (1948) . The eddy viscosity for
shear layers in a stratified flow is defined by
Evo - « 2 Ml - jf> <3'6'
where K = von Karman's constant («0.4)
z = distance above the bottom
U* = shear velocity
h = depth of flow.
3.5
-------
To account for stratification
E..
Evo
(1 + 10 Ri)'1/2 (3.7)
and
£*-- (1 + 3.33 Ri)-3/2 (3-8)
So
where E and K are the values of eddy viscosity and diffusivity
for the neutral case. This model was found to work reasonably well
for most flows, but problems occurred in situations in which mixing
through and below the thermocline occurred. A site-specific method
takes the form
p- = (1 + b Ri)"n for Ri < 1 (3.9)
vo
and
TT-= (1 + b)'n for Ri > 1 (3.10)
vo
where b and n are empirical coefficients.
Wang (1982) used a form of the Munk and Anderson
formula for modeling flows around a circular island with a shelf
along the shore. Two models were run. In the first run, the shelf
sloped from 40 to 200 m at the shelf break, and in the second model
the depth was held constant at 100 m. The eddy viscosity and
diffusivity were modeled by
3.6
-------
Ey = 5 + 50 (1 + 10 Ri)"1/2 (3.11)
Kv = 50 (1 + 3.33 Ri)"3/2 (3.12)
respectively. If one compares the results for a case in which Ri
— * 0, the values obtained compare with the values for the surface
layer in Table 3-1. This is to be expected in the shallow depths
being modeled.
Blumberg (1977) , in a paper on estuarine circulation, used
another formulation for the definition of the eddy diffusivity:
where k- = empirical constant («0.10 for tidally averaged flow)
z = distance above the bottom
h = depth of flow.
For diffusion adjusted for stratification:
1/2
v i/2 ,2 n z. 2idUi , , Ri \ (3
Kv = kl z (1 " h> ldzl ( l RiJ (
and
E.. = Kt, (1 +
Ri is the critical Richardson number beyond which turbulence is
c
suppressed by stratification. There is some disagreement as to
what this value should be. Blumberg (1977) recommends a Ric value
of 10, which seems high when compared with other recommendations.
Geyer (1973) did extensive studies in the Fraser estuary and found
3.7
-------
that a Ri of 0.25 was the point beyond which mixing is suppressed.
c
It is possible that internal waves cause mixing to occur at higher
Richardson numbers, but the maximum observed value was
approximately 1. A value of 0.25 agrees with other experimental
values of Ri and is the preferred choice for a critical Richardson
c
number.
K will approach zero as one nears the surface (Blumberg
1977). This does not agree with the values of Table 3.1, which
show the value of KV approaching 50 cm2/s. This large eddy
diffusivity is probably caused by wind stress on the surface and
will vary with the magnitude of the stress. It is possible to
avoid the eddy diffusivity approaching zero at the surface by
measuring z from the surface in the upper regions.
A high eddy diffusivity is found in the upper region, where
surface stresses cause active mixing. The diffusivity decreases
as one reaches the middle region and remains relatively constant.
Near the bottom the diffusivity again increases because of the
bottom stress, although it is not as great as near the surface.
Zero-equation models are simple and have been used
extensively in engineering, but the method is rather crude.
Although history and transport are not taken into account, quite
good results have been obtained in some cases. If a reasonable
relationship can be found for the interior region, this method is
the simplest to use.
3.2 ONE-EQUATION MODELS
One-equation models require that one additional partial
differential equation be solved to obtain a closure for the
turbulence equations. The eddy viscosity is set proportional to
a velocity scale and a length scale, and a partial differential
3.8
-------
equation is developed for the velocity scale. The length scale is
determined by empirical means. The use of the differential
equation allows the model to include the history and transport
effects of the turbulence. The most common velocity scale is
turbulent kinetic energy (k). The exact turbulent kinetic energy
(TKE) equation can be obtained through manipulation of the
Navier-Stokes equations (Rodi 1984)
u.u.
au. _ au.au.
where P = pressure
p = density
/3 = volumetric expansion coefficient
gi = gravity.
The exact equation involves the turbulence
correlations and requires that several model assumptions be made
to solve the equation. The diffusion flux is set proportional to
the gradient of k
u
u
•i v 2 T p' o. ax.
ni I
The dissipation is modeled by
c
CD L
3.9
-------
Using these assumptions and the eddy viscosity/diffusivity
expressions for U| uj and u( 0, the equation can be written as
r jal I jal I
ak . „ dk _ _a_ ,S/. ak_, F ,_i
at + ui ax~ - aXi (ak ax/ + Ev (d*.
ax.
where a, and CD are empirical constants, and a^ is the turbulent
Prandtl or Schmidt number.
Rodi (1984) uses the equation
Ev = CE Vk L (3.18)
to define the eddy viscosity, where C is an empirical constant,
and L is a length scale. It has been found that C_ and C_ « 0.08,
Ci U
and a. « 1. Once the eddy viscosity has been determined, the eddy
JC
diffusivity can be calculated using the Prandtl or Schmidt number,
and the mean flow equations can be solved.
A one-equation model includes transport of turbulence
and is therefore preferred over a zero-equation model in areas
where transport is important. The difficulty with this method li^s
in the fact that a turbulence length scale must still be
determined. In regions where a length scale can be clearly
defined, the one-equation model can work quite well. A boundary
layer is an example. In areas where no length scale can be clearly
defined, the method tends to break down, and two-equation models
would be a better choice.
Several methods have been developed to determine the length
scale. An empirical formula that uses one length scale when
3.10
-------
turbulence production dominates and a second scale when convection
is dominant is described by Rodi (1984). Bernard's method seems
to work well, but in some situations different constants are
required to accurately describe the flow field. A common method
to determine the length scale, as originally proposed by von
Karman, uses local derivatives of the velocity- Gawain and
Pritchett (1970) developed different methods for the length scale,
but there are no data currently available to validate the models
over a wide range of flows. A different one-equation model was
developed by Bradshaw (1978). In this model a transport equation
for the shear stress, uv, was solved instead of using the eddy
viscosity concept.
All these methods are fairly complicated, and some require
considerable computer time because they are iterative methods.
Also, all these methods were developed for wall-boundary layer
flows and would require modification if they were used in a free-
flow situation. The simplest methods are those based on the ideas
of von Karman, which should be considered before attempting the
more complicated methods of Bernard or Bradshaw (1978), in which
no length scale is clearly defined.
3.3 TWO-EQUATION MODELS
These models require the solution of two partial different
equations to obtain a closed solution of the mean velocity field
equations. In the previous section, the TKE equation was presented
as a velocity scale; here, an equation will be defined for the
length scale. This allows the use of transport and history effects
in the length scale as previously used in the velocity scale.
One of the most common ways to define the length scale is to
use the rate of TKE dissipation (e), resulting in the k-e model.
3.11
-------
Others disagree with the term e and have developed a method using
the term kL where k is the TKE, and L is the length scale. These
two methods are discussed below.
The transport equation for the rate of TKE dissipation is
f * «i - (f £ + cu f - C2c r <3
where a , C. , C_ , and C- are empirical constants, and P and G
are turbulence production due to stress and buoyancy, respectively.
The constants C. coe' an(* ae nave ^een found to be 1.44, 1.92,
and 1.3, respectively. The value for C_ varies depending on the
flow field; when G is a source term, C. — 1.0, and when G is a
sink term, C_ «. 0. There are some questions as to the validity
of the constants, but they have proved accurate in many test
situations and are a good starting point.
After the determination of k and e, the eddy viscosity and
diffusivity can be found from
E - c <3-20'
Kv - L (3.2!)
where GU is an empirical parameter assumed constant (A 0.09), and
o. is the Prandtl or Schmidt number.
There are several problems with the e equation, which is
based on the idea of isotropic dissipation. While this may be true
sometimes, it does not allow for variation of the eddy viscosity
3.12
-------
and diffusivity in different directions, and it is known that the
horizontal values of the eddy viscosity/diffusivity are several
orders of magnitude greater than the vertical components. The
isotropic restriction has been eased by some researchers who have
varied the eddy viscosity/diffusivity in the vertical and
horizontal directions. Although this does seem to work, it does
not accord with the strict definition of the equations.
Mellor and Yamada (1982) disagreed with the use of the energy
dissipation rate as the length scale and have developed an equation
using kL:
akL M <9kL d , ff - g- .
aT + ui aT = d^7 (^ s^ L
L (P + G) -
where S£/ E^ E^f and B^ = 0.2/ 1.8, 1.33, and 16.6, respectively,
and K is von Karman's constant («0.4). Rodi (1987) states that the
argument on the relative merits of either equation is rather
academic, because both equations are fairly empirical, and with the
correct constants these methods should perform satisfactorily.
Mellor and Yamada (1982) believe that using e is wrong because it
describes small-scale turbulence, and the length scale should be
based on large-scale turbulence. Either method should work
satisfactorily if the correct constants are used.
A problem arises when the effects of stratification must be
taken into account. It has been found in some cases that the value
of C in Equation (3.20) is not really constant but varies with the
buoyancy effects. A model known as the extended k-e model
(algebraic stress/flux model) is given in Rodi (1987) as a method
3.13
-------
to solve the stratification problem
(3.23)
1 k
10
(3.24)
where the empirical constants, C. , C_, C.., C_,, and R are 1.8,
L £ L
-------
4. MODEL DESCRIPTION
4.1 TEMPEST/FLESCOT MODEL
The TEMPEST model is an unsteady, three-dimensional finite
difference model (Trent et al. 1989; Onishi et al. 1985a,b) cast
in cylindrical or Cartesian coordinates. It simulates flow,
turbulent kinetic energy (TKE), heat transfer, and general
dissolved mass transfer (e.g., salt). The model has been con-
tinuously updated and modified to expand its applicability to a
wide range of hydrodynamic problems. For example, it is currently
being modified to use an orthogonal curvilinear coordinate system.
It can compute local (heterogeneous) isotropic eddy viscosity and
dispersion coefficients by using the k-e TKE and dissipation model
if this option is chosen. Otherwise, it uses eddy viscosity as
input data.
The model can calculate these variables with either the
dynamic pressure approach or the hydrostatic pressure assumption
based on
• conservation of mass for fluid (the continuity equation)
• conservation of momentum (the Navier-Stokes equations)
• conservation of TKE (k-e model)
• conservation of heat (the first law of thermodynamics)
• conservation of mass for constituents including salt.
The FLESCOT model is a sediment-contaminant transport
version of the TEMPEST model. Thus, it is also an unsteady, three-
dimensional finite difference model (Onishi and Trent 1982).
FLESCOT simulates time-varying movements of flow, TKE, salt, heat,
sediment, and dissolved and sediment-sorbed contaminants (e.g.,
radionuclides, heavy metals, pesticides, PCBs, etc.) in surface
4.1
-------
waters. In addition to the equations described for the TEMPEST
model, FLESCOT has equations describing sediment-contaminant
transport and accumulation on the bottom. Sediment movement and
particulate contaminant transport in FLESCOT are modeled
separately for three sediment size fractions or sediment types.
The sediment transport submodel includes the mechanisms of
advection and dispersion of sediments, fall velocity and cohesive-
ness, deposition on the seabed or ocean bottom, erosion from the
bed (bed erosion and armoring) , and sediment contributions from
point/nonpoint sources and subsequent mixing. This submodel also
calculates changes in bed conditions, including bed elevation
changes caused by scouring and/or disposition, and gives a three-
dimensional distribution of sediment sizes within the bed.
FLESCOT also includes wave mechanisms to enhance bottom shear
stress and increase suspended sediments in shallow water, based on
Grant's model (Grant and Madsen 1979).
Dissolved contaminants interact both with sediments
in motion (suspended and bed-load sediments) and with stationary
sediments on the seabed or ocean bottom. To account for the
interactions, the transport of sediment-attached contaminants in
FLESCOT is solved separately for each sediment size fraction. The
model includes the mechanisms of advection and dispersion of
particulate contaminants, adsorption/desorption of dissolved
contaminants with sediment, chemical and biological degradation (or
radionuclide decay, if applicable) of contaminants, deposition of
particulate contaminants on the bed or erosion from the bed, and
contaminantion from point/nonpoint sources and subsequent mixing.
The three-dimensional distributions of the particulate contaminants
within the bed are also computed.
For simplicity, we use Cartesian coordinates to express
governing equations of the model.
4.2
-------
4.2 HYDRODYNAMIC EQUATIONS
Continuity equation
dy. , dv , dw _ n /AM
ax + ay + dz ~ ° t4-1)
where u, v, and w = velocity components in the x-, y-, and
z- directions, respectively; and,
x, y, and z = longitudinal, lateral, and vertical (upward)
directions, respectively.
Momentum equations
For the x-component,
_ 1 ap_ xa_ , au,
r = + e
_ _
at ax ay <9z x = P ax
For the y-component,
d\/_ auv avv avw ,f_ I5P,^_/t5v,,a_
at dx dy dz Ty p dy dx iex axj ay
For the z-component,
aw auw avw aww ,f= lap n , 5_ / , 5w >
at ax ay az z p az y ax ^ex ax;
4.3
-------
where f = 2n (w cos 0 v sin 0) (Coriolis force in x-direction)
A
f = 2n u sin 0 (Coriolis force in y-direction)
f = 2n u cos 0 (Coriolis force in z-direction)
F , F = bottom shear stress and surface (wind) shear stress,
A A
respectively, in x-direction
bottom shear stress and surfc
respectively, in y-direction
F , F s = bottom shear stress and surface (wind) shear stress,
F = shear stress along the vertical solid boundaries
g = acceleration due to gravity
P = pressure
p = water density
0 = geographical latitude
n = rate of earth's rotation (= 7.29 x 10"5 s l)
V V ez = eddy viscosity components in the x-, y-, and z- directions,
respectively.
Note that the hydrostatic pressure assumption is not made and
that the vertical momentum equation retains all terms of the flow
acceleration, Coriolis, eddy viscosity, and frictipn, in addition
to pressure and gravitational terms.
For wind shear stresses,
= c* p W2 sin
(4.5)
FyS= c* pa Wa2 cos * (4.6)
*
where c = wind stress drag coefficient
Wa = wind speed measured at some specified height (such as 10 m)
p. = air density
a
* = angle between the wind and current directions.
4.4
-------
Similarly, bottom friction is commonly expressed by
c b / 2 2
FX = apu -J u + v (4.7)
Fy = apv J u2 + v2 (4.8)
FZ = apw J u + w or apw J v2 + w2 (4.9)
where coefficient "a" is a dimensionless drag coefficient that can
be correlated to Manning and Chezy coefficients. If the Chezy
coefficient, C, is used, the above expressions become
(4.10)
(4.1D
Eddy viscosity in the Navier-Stokes equations must be either
assigned as input data or computed internally, based on the k-e
model .
In addition, the following equation of state relates fluid
density to constituent concentrations (e.g., water temperature,
salinity, and possibly sediment concentrations) :
P = R (T, S) • (4.12)
where R = functional relationship
T = water temperature
S = salinity
4.5
-------
Note that the water density is not a function of pressure when the
incompressibility assumption is made.
The following kinematic boundary condition on the water
surface is used to express a water surface
K = u <%L v K , w (4
dt u <9x v dy + w I*
where f = water surface elevation.
TEMPEST/ FLESCOT has built in two options of the pressure
calculations; one is to assume the commonly used hydrostatic
pressure assumption, and the other is not to assume the hydrostatic
pressure assumption, thus solving Equation (4.4) with all the terms
included. Solving the continuity and momentum equations is
difficult, especially the iterative pressure calculation in
Equations (4.2), (4.3), and (4.4). This process (dynamic pressure
calculation approach) of not assuming the hydrostatic pressure
assumption is used by a few hydrodynamic models, such as
TEMPEST/FLESCOT (Trent et al. 1989; Onishi and Trent 1982), FLOWER
(Eraslan et al. 1983), and REMIXCS (Freitas et al. 1985).
Fortunately, most flows occurring in the natural environment
have very small vertical velocities and essentially negligible
vertical acceleration, as compared to the gravitational term.
ThuSj the hydrostatic pressure assumption can be made in most
cases. The adaptation of the hydrostatic pressure assumption
reduces the vertical equation of motion, Equation (4.4), to the
following simple form:
(4.14)
4.6
-------
where the assumption is made that the temporal and spatial vertical
velocity accelerations and shear stress changes (eddy viscosity
terms) are much smaller than the gravitational acceleration. Thus,
under the hydrostatic pressure assumption, TEMPEST/FLESCOT uses
Equation (4.14) instead of Equation (4.4). The vertical velocity
in the Coriolis force may also be dropped in this case. Under the
hydrostatic pressure assumption, the kinematic boundary condition
Equation (4.13) is further simplified by dropping the water surface
gradient terms in the x- and y- directions, resulting in the
following equation:
f«w (4.15)
4.3 HEAT TRANSFER
Heat transfer equation
fe
where C = specific heat
qT = rate of heat generation or dissipation including heat loss or
gain through water surface and/or ocean bottom
p = fluid density
a , a , a = thermal conductivity in the x, y, and z directions,
x y z
respectively.
4.7
-------
4.4 MASS TRANSPORT
Dissolved mass transfer ecfuation
dC , dC dC dC _ d_ /, dC\
<3y y <3y' dy v z dz' Mc v
where C = concentration of mass (e.g., salt)
k , k , k = dispersion coefficients in the x, y, and z directions,
x y z
respectively
q = rate of mass generation or dissipation.
These equations or corresponding equations in other coordinate
systems are solved by many three-dimensional computer codes to
obtain distributions of velocity, water temperature, salinity, and
other constituents.
4.5 TURBULENCE EQUATIONS
Modeled transport equations for TKE, k, and dissipation of
TKE, e, are solved to obtain the turbulent effective viscosity, ju .
As discussed in Chapter 3, the modeled forms of the equations are
Turbulent kinetic energy
CTk
where £. = ^
Sk =Pk+Gl
4.8
-------
Shear production
, / C7U OV \ C t OU C/W \ L / OU CW i t ^ / A 1 Q \
ay ax' ^az alT' ^az ay*' ' \ • i
Buoyant production
Dissipation of turbulent kinetic energy
If
/c \ , _ /p \ , _ c
ax Ke ax; ay ^e ay; az ^e az;
F (Se pCe2e)€
where
= C£lPk
where M = dynamic viscosity
4.9
-------
The turbulent viscosity, MT/ is computed using the Prandtl-
Kolmolgorov hypothesis:
^ = CM p k*/e (4-22)
Recommended turbulent model constants (Jones and Launder 1973) are
a = 0.9 C . = 1.44 C = 0.09
T c 1 M
CT, = 1.0 C_=1.92
k £2
"e - 1-3 Ce3 = X-44
4.6 SEDIMENT TRANSPORT EQUATION
The governing equation for the transport of the j-th sediment
in the FLESCOT model is
^j d_ . „ , d_ , p , 6_ r, )c -,
<3t <3x j <9y j (9z sj j
a oU • - C7L • , C?L • Ri D1
= ax ^ x dx ' + dy~ ^ y ay ^ + 5z ^ z az ^ + ^ h " h ' + qcj ^ ' '
where Cj = concentration of j-th sediment
h = water depth
qcj = rate of sediment generation
SDJ = j-th sediment deposition rate per unit bottom surface
area
SRJ = j-th sediment erosion rate per unit bottom surface area
W8, = fall velocity of the j-th sediment.
4.10
-------
To obtain the sediment concentration, sediment erosion and
deposition rates, calculate SRj and SDJ for each sediment size
fraction separately. For noncohesive sediment (e.g., sand), if the
amount actually being transported is less than the flow can carry
for given hydrodynamic conditions, the flow will scour noncohesive
sediment from the seabed and ocean bottom. This process will
increase the sediment transport rate with the rate based on the
difference between the sediment transport capacity of the flow and
the actual sediment transport rate. The process occurs until the
actual sediment transport rate becomes equal to the carrying
capacity of the flow or until all the available bottom noncohesive
sediment is scoured, whichever occurs first. Conversely, non-
cohesive sediment is deposited with the deposition rate again based
on the difference between the actual sediment transport rate and
sediment transport capacity of the flow. Because of the simplicity
of the formulation, FLESCOT currently uses DuBoy's formula (Vanoni,
1975) to estimate the noncohesive sediment transport capacity of
flow. Partheniades' (1962) and Krone's (1962) formulas are used
to calculate cohesive sediment erosion and deposition rates.
FLESCOT divides the seabed and ocean bottom into a number of
bottom layers. Contaminant distributions associated with each
sediment size fraction within each bottom layer are also obtained
by keeping track of the amount of contaminants removed from or
added to each bottom layer during the simulation period. The
mechanisms of erosion and deposition of contaminated or clean
sediment, and direct absorption/desorption occurring between dis-
solved contaminant and bottom sediment are used in determining
contaminant distributions in the bottom.
4.11
-------
4.7 DISSOLVED CONTAMINANT TRANSPORT EQUATION
A governing equation for a dissolved contaminant is
<9G A a A
|r (vGJ + f- (wGJ
dt dx v w' dy v w' dz v w'
. (9G ,, <3G
ax ' dy y dy dz z dz
i (4-24)
IzTj(l-POR)DjKBj (KdjGw GBj)
J
h Tj (l-POR) Dj KBj (Kdj Gw GBj,
J
where Dj = diameter of j-th sediment size fraction
KBJ, K^j = transfer rate of contaminants for adsorption and
desorption, respectively, with j-th nonmoving sediment in bottom
Kdj/ Kd'j = distribution (or partition) coefficient between
dissolved contaminant and particulate contaminant associated with
j-th sediment for adsorption and desorption, respectively
GBj = particulate-contaminant concentration per unit weight
of sediment in j-th sediment size fraction in the bpttom
GJ = particulate-contaminant concentration associated with
j-th sediment (radionuclide activity or weight of contaminant) per
unit volume of water
Gw = dissolved-contaminant concentration (radionuclide
activity or weight of contaminant) per unit volume. of water
Gwo = constant concentration of dissolved contaminant
FOR = porosity of bottom sediment
Qw = lateral influx or other source strength of dissolved
contaminant
4.12
-------
7j = specific weight of j-th sediment
7 = radionuclide decay, or chemical and biological degradation rates
of contaminant, where applicable.
The particulate-contaminant concentration associated with sediment
in the water column, Gj, in Equation (4.24) is expressed in terms
of the contaminant weight or radionuclide activity per unit volume
of water, instead of per unit weight of sediment. The printout of
particulate radionuclide concentration, however, is converted to
contaminant weight or radionuclide activity per unit weight of
sediment.
4.8 SEDIMENT-SORBED CONTAMINANT TRANSPORT EQUATION
A governing equation of a sediment-sorbed contaminant
associated with j-th sediment is
A
IT
-------
(4.12) is first used to update the water density, q, with known
values of T and C. Then Equations (4.2) through (4.4) are used to
calculate velocity components, u, v, and w. Equation (4.1),
together with some residual terms in Equations (4.2) through (4.4),
are then used iteratively to update u, v, w, and P to satisfy the
fluid mass balance. This iteration is needed because the values
of u, v, and w calculated by Equations (4.2) through (4.4) use
values of P and q at the previous time step. These velocities do
not necessarily satisfy fluid mass continuity at the present time
step with a current pressure value. Corrected velocity values,
through the fluid mass balance iteration procedure, are then used
to calculate T and C with Equations (4.16) through (4.25). This
procedure is repeated at every time step.
Using the hydrostatic pressure assumption, the
TEMPEST/FLESCOT model solves the unknowns faster and more easily.
As before, q is first updated using the known values of T and C.
Then by using the horizontal momentum equations [Equations (4.2)
and (4.3)], the horizontal velocity components, u and v, are calcu-
lated with known values of P and q. By using calculated u and v,
the vertical velocity, w, is then computed using the continuity
Equation (4.1) alone. The new water surface is then computed by
Equation (4.15). A new pressure distribution is calculated from
Equation (4.14) with the new water surface elevation. The
calculated u, v, and w are then used to calculate T and C with
Equations (4.16) through (4.25).
This much simpler hydrostatic pressure approach to calculate
velocity and pressure is much faster (one to two orders) than using
the dynamic pressure approach. The vertical velocity calculation
is decoupled from the horizontal velocity calculations, and there
is no need to iterate to obtain the pressure value.
However, the hydrostatic pressure assumption is less
4.14
-------
applicable to an area where vertical velocity and acceleration are
important. Potential computational errors with the hydrostatic
pressure assumption are potentially greater than those associated
with the dynamic pressure assumption because any computational mass
balance errors in u and v calculations with Equations (4.2) and
(4.3) will be absorbed by the vertical velocity, water surface, and
pressure calculations without corrective (iterative) actions taken
by the dynamic pressure approach.
4.15
-------
5. MODEL APPLICATIONS
The TEMPEST/FLESCOT code was applied to a two-dimensional,
3000 in-deep, rectangular basin to determine its applicability to
a deep-ocean environment. It was first applied to calculate
vertical eddy viscosity; then to calculate flow, water temperature,
sediments, and radionuclide distributions.
5.1 EDDY VISCOSITY DISTRIBUTIONS
The code was applied to a 12,000-m-long, 3000-m-deep
rectangular basin under three stratified and nonstratified
conditions. For 3 sets of eddy viscosity calculations, we used a
dynamic pressure approach, even though a hydrostatic pressure
approach can also be used. Note that we used computational cells
near the bottom with very small vertical increments (0.1 m) to
obtain detailed distributions of near-bottom flow, turbulence, and
eddy viscosity. Because longitudinal grid spacing is fixed at 500
m, there is a large aspect ratio of computational cells near the
bottoms. Their impact on the k-e model is not clarified in this
study. For each simulation case, a constant of 0.1 m/s was used
for initial and boundary velocity distribution values.
5.1.1 Nonstratified Case (Case 1)
The first case was a nonstratified condition with water
temperature of 10 °C and salinity of 32 °/oo over the entire flow
region. Figures 5-1 to 5-3 show distributions of calculated
turbulent kinetic energy (TKE), turbulent energy dissipation, and
eddy viscosity for 2-days simulation over the entire flow region.
Contour levels for TKE (Figure 5-1) are 1 x 10"10, 1 x 10"8,
1 x 10"6, 5 x 10"6, 1 x 10"5, 1.5 x 10"5, 2 x 10'5, and 3 x 10"5 m2/s2.
Contour levels for turbulent energy dissipation (Figure 5-2) are
I x 10"15, 1 x 10"12, 1 x 10"10, 1 x 10"9, 5 x 10"9,1 x 10"8, 5 x 10"8,
5.1
-------
1 x 10"7, 5 x 10"7, and 1 x 10"6 m2/s3. For eddy viscosity, the
contour levels used in Figure 5-3 are 1 x 10"4, 5 x 10"4, 1 x 10"3,
5 x 10"3, 1 x 10"2, 5 x 10'2, 1 x 10" V 5 x 10~1, 1, 5, 10, 15, and 20
Pascal-seconds. Note that the values are highest near the bottom
and rapidly decrease as the distance from the bottom increases.
Vertical ctontour lines (Figures 5-1 and 5-3) are due to
extremely small values assigned at the upstream boundary. Figures
5-4 to 5-6 show these values within approximately 70m of the
bottom. Contour levels for TKE (Figure 5-4) and turbulent energy
dissipation (Figure 5r5) are the same as used in Figures 5-1 and
5-2, respectively. Contour levels for the eddy viscosity
(Figure 5-7) are the same as those in Figure 5-3, except that the
levels of 1 x 10"4 and 5 x 1Q"4 Pascal-seconds were eliminated.
Figures 5-7 to 5-9 show vertical distribution of TKE, turbulent
energy dissipation, and eddy viscosity at 10,000 m downstream.
Figures 5-4, 5-5, 5-7 and 5-8 reveal that the TKE is generated most
at the bottom and that the turbulent energy is dissipated most at
the bottom. Figures,5-6 and 5-9 indicate that the eddy viscosity
calculated by the k-e model increases originally with the vertical
distance from the bottom, almost linearly, until it peaks out about
30 m from the bottom. Values of eddy viscosity sharply dropped to
levels near that of molecular viscosity about 65 m from the bottom.
The TKE and eddy viscosity reduce their calculated values
sharply with vertical distance beyond approximately 45 m from the
bottom. This is because the vertical gradient of the calculated
longitudinal velocity above 45 m is almost zero (the longitudinal
velocity above 45 m reaches its uniform value of 0.1 m/s), thus no
Shear production of TKE due to the velocity gradient (see Equation
4.19) occurs above 45 m, resulting in sharp reduction of the TKE
and eddy viscosity values above that height. The maximum
calculated eddy viscosity after a 2-day simulation is approximately
20 pascal-seconds (or 200 cm2/s) occurring about 30-45 m above the
5.2
-------
bottom. This value is higher than those shown in Table 3-1 and may
be due to the lack of stratification which suppresses the TKE, and
thus eddy viscosity.
As discussed in Case 3 with stratification near the bottom,
the calculated eddy viscosity values near the bottom are very close
to the value of 1 cnf/s reported in Table 3-1. However, there are
no specific values available for comparison. Rresulting calculated
velocity distribution (Figure 5-10) appears as well-developed,
flattened velocity profiles.
5.1.2 Stratified Cases (Cases 2 and 3)
We examined two stratified cases where stratification was
imposed by water temperature distribution, not by salinity
distribution, both of which are shown in Tables 5-1 and 5-2.
TABLE 5.1. Initial Water Temperature and
Salinity Distribution for Case 2
Distance Water
from Bottom (m) Temperature ( C) Salinity f°/oo)
2960 - 3000 16 32
2930 - 2960 13 32
150 - 2930 10 32
120 - 150 6 32
0 - 120 2 32
5.3
-------
TABLE 5.2. Initial Water Temperature and
Salinity Distribution for Case 3
Distance Water
from Bottom (m) Temperature f°C) Salinity f°/oo)
2960 - 3000 16 32
2950 - 2960 13 32
7.5 - 2950 10 32
6.0 - 7.5 6 32
0-6 2 32
Calculated TKE, turbulent energy dissipation, and
corresponding eddy viscosity over the entire flow field for Case 2
are shown in Figures 5-11 to 5-13. Those within 70 m from the
bottom are shown in Figures 5-14 to 5-16. Contour levels for
Figures 5-11 to 5-16 are the same as those used in corresponding
Figures 5-1 to 5-6 for Case 1. Vertical distributions of TKE,
turbulent energy dissipation and eddy viscosity at 10,000 m
downstream are shown in Figures 5-17 to 5-19. These figures show
the small reduction of these values from those obtained in
nonstratified Case 1), because of the stratification. For Case 2,
the maximum calculated eddy viscosity is approximately 17
2
Pascal-seconds (170 cm /s) occurring at 45 m above the bottom.
This is due to the existence of the density gradient which reduces
both the TKE and the dissipation of turbulent kinetic energy
through negative buoyancy production (see Equations 4.18, 4.20, and
4.21). The net result also reduces eddy viscosity. However, since
the TKE is generated near the bottom and the initial water
temperature was constant (2°C) at about 200 m from bottom, effects
of stratification on TKE, turbulent energy dissipation, and eddy
viscosity were very small in comparison to Case 1. Because of a
density gradient about 120-150 m from bottom, calculated
5.4
-------
longitudinal velocity did not become uniform until approximately
250 m above the bottom. This results in more gradual variations
of TKE and eddy viscosity in the bottom 78 m (Figures 5-17 and 5-
19) , as compared to Case 1 (Figures 5-7 and 5-9). The calculated
velocity distribution is shown in Figure 5-20.
To further test stratification effects on TKE dissipation,
and thus on eddy viscosity, PNL applied the code to Case 3
conditions having water temperature gradients - including those at
6 to 7.5 m above the bottom (see Table 5-2). Calculated
distributions are shown in Figures 5-21 to 5-30. Contour levels
of TKE used in Figures 5-21 and 5-24 are 1 x 10'10, 1 x 10"9,
1 x 10*, 5 x 10"8, 1 x 10"7, 5 x 10"7, 1 x 10"*, and 5 x 10"6 m2/s2;
Contour levels of turbulent energy dissipation in Figures 5-22 and
5-25 are 1 x 10"13, 1 x 10"12, 1 x 10'11, 1 x 10'9, 5 x 10'9, 1 x 10"8,
5 x 10"8, and 1 x 10"7 m2/s3. Contour levels of the eddy viscosity
(Figures 5-23 and 5-26) are 1 x 10"4, 5 x 10"4, 1 x 10"3, 5 x 10"3,
1 x 10"2, 5 x 10"2, 1 x 10"1, and 2 x 10"1 Pascal-seconds. Since Case
3 has a density gradient (6 to 7.5 m above the bottom), where the
activities of the TKE and turbulent energy dissipation are the
greatest, this case clearly demonstrates the suppression effect of
t
stratification (buoyant production term, Equation 4.20) on TKE,
turbulent energy dissipation, and eddy viscosity. The reduction
of TKE and eddy viscosity are particularly significant. Note the
significant reduction of computed eddy viscosity in this case,
having the maximum eddy viscosity values of approximately
0.2 Pascal-second (2 cm2/s) at 1 m above the bottom, as compared
with the maximum value of 20 Pascal-seconds (200 cm2/s) for the
nonstratified Case 1. The calculated eddy viscosity of 2 cm2/s for
this case agrees well with the 1.0 cm2/s vertical eddy diffusivity
estimated by the U.S. AEC (1975) at 2000*m deep (see Table 3-1).
The density gradient near the bottom also changed the computed
velocity distribution near the bottom, as follows: the
5.5
-------
longitudinal velocity in this case reaches a uniform value of
approximately 0.2 m/s closer to (about 20 m above) the bottom,
Within this bottom 20 m, however, exists a velocity variation
(especially within the bottom 0.5 m) that is much more gradual with
vertical distance, displaying a more parabolic character than the
flattened turbulent velocity distribution seen in Cases 1 and 2.
v
The model demonstrates that the density gradient must be
considered in calculating eddy viscosity, and that the k-e model
reproduces, at least qualitatively, the vertical eddy viscosity
variation in a deep-ocean environment. Although we have yet to
apply the k-e model to a three-dimensional deep water case, the
TEMPEST code has been applied successfully to three-dimensional
thermal plume cases to calculate turbulence with the k-e model.
5.2 RADIONUCLIDE DISTRIBUTIONS
5.2.1 Simulation Conditions
Using the eddy viscosity calculated in Cases l and 2, three
additional cases (Cases 4 through 6) were simulated to obtain the
following distributions:
velocity
water temperature
suspended sand concentration
suspended silt concentration
suspended clay concentration
dissolved radionuclides
radionuclide sorbed by sand
radionuclide sorbed by silt
radionuclide soibed by clay.
Note that the code also calculated the change on bottom conditions;
that is, the radionuclide distribution sorbed by the three
fractions (sand, silt, and clay) of bottom sediment.
In all three cases (Cases 4 through 6), it was assumed that
5.6
-------
there were no radionuclides in the water column initially, and the
radionuclides originally existed only in bottom sediments of the
500-m reach between 2000 and 2500 m from the upstream end of the
study area. Conditions imposed in Cases 4 through 6 are shown in
Table 5-3. Note that to save computation time, the hydrostatic
pressure approach was selected for Cases 4 through 6.
Eddy viscosity at each computational cell for each case was
assigned, as follows, to allow for a nonisotropic field:
Vertical = molecular viscosity + eddy viscosity, as
calculated in Case 1 or 2
Longitudinal = molecular viscosity + eddy viscosity,
of Case 1 or 2, + oceanic lateral or
longitudinal eddy viscosity (assumed to
be 1000 Pascal-seconds).
Nonisotropic dispersion coefficients were estimated by the
corresponding eddy viscosity and Schmidt/Prandtl numbers (which are
assumed to be 0.71) by using Equation 3.4.
5.2.2 Simulation results of Stratified Case 4
As stated above, by using calculated eddy viscosity, the
code then simulated radionuclide transport phenomena. Because of
the very small vertical grid spacing near the 3000 m bottom-depth,
we used a 0.1-second time-step for Cases 4, 5 and 6 to ensure
computational stability. The calculated results for radionuclide
concentrations are those obtained after a 6-hr simulation.
The model predicted that a small portion of radionuclides,
sorbed by bottom sediment, leached out to the overlying water
volume. Once in the water-column, some of the leached, dissolved
nuclides were then transported downstream by current (s) . Moreover,
some were absorbed by suspended silt and clay. The model also
predicted that the ratio of the radionuclides sorbed by suspended
silt and clay is equal to the ratio of distribution coefficients
5.7
-------
TABLE 5.3. Simulation Conditions and
Model Parameters, Cases 4 through 6
Case 4 Case 5 Case 6
Stratified Condition Same as Same as Same as
Case 2 Case 2 Case 1
Sediment Diameter (mm)
Sand 0.25 0.25 0.25
Silt 0.016 0.016 0.016
Clay 0.002 0.002 0.002
Specific Weight (g/cm3)
Sand 2.65 2.65 2.65
Silt 2.65 2.65 2.65
Clay 2.65 2.65 2.65
Settling Velocity (m/s)
Sand 5 x 10~l 5 x lO* 5 x 10^
Silt 5 X lO* 5 X 10* 5 X lO*
Clay 1 x 10"6 1 x 10"6 1 x 10"6
Initial Concentration
in water-column (mg/L)
Sand 000
Silt 0.1 100 100
Clay 0.1 100 100
Initial Fraction
in bottom
Sand 0.2 0.2 0.2
Silt 0.6 0.6 0.6
Clay 0.2 0.2 0.2
Radionuclide Concentration
at bed source (mCi/g)
Sorbed by sand 0.001- 0.001 0.001
Sorbed by silt 111
Sorbed by clay 555
Distribution Coefficient,
Kd, (ciir/g) with
Sand 0.02 0.02 0.02
Silt 20 20 20
Clay 100 100 100
Adsorption/Desorption
Rate, Kj, (1/s^
Sand 1 x lol 1 x 10JJ 1 x 10"!
Silt 1 x lot 1 x lot 1 x'10":
Clay 1 x 10"* 1 x 10"6 1 x 10'6
5.8
-------
associated with silt and clay, and that the originally "clean"
bottom sediments, downstream of the contaminated bottom portion,
then adsorb a fraction of the radionuclides in the water-column
back into the bed.
Figures 5-31 to 5-33 show the predicted distribution for
dissolved, suspended-silt-sorbed and suspended-clay-sorbed
radionuclides in the water column up to 40 m above the bottom.
Note that the unit of dissolved concentration is pCi/L, while the
sediment-sorbed radionuclide concentrations are in pCi/g.
Table 5-4 shows the predicted radionuclide distribution in
the top 3-cm bed layer in pCi/g. Note that the bottom sediment
between 2000 to 2500 m is the original radionuclide source.
TABLE 5.4. Calculated Radionuclide Concentration
Sorbed by Sediment in the Top 3-cm Bed Layer After 6 h Simulation
Downstream
Distance
m
Sand
Radionuclide Concentration. pCi/g
Composite
,Sediment
Silt
Clav
2000-2500
(source)
2500-3000
3000-3500
3500-4000
4000-4500
4500-5000
9.993x10-4
5.99x10-13
2.25x10-13
9.98X10-14
4.47x10-14
1.96x10-14
9.9998x10-1
1.60X10-11
5.99x10-12
2.64x10-12
1.19x10-12
5.24X10-13
4.99996
3.00x10-11
1.12x10-11
4.94x10-12
2.23x10-12
4.82x10-13
1.60018
1.59x10-11
5.89x10-12
2.59x10-12
1.17x10-12
5.14x10-13
5.2.3 Simulation Results of Stratified Case 5
*
The only difference between Cases 4 and 5 is that the suspended
sediment concentration for Case 5 has values that are 1000 times
5.9
-------
higher than in Case 4. Although a sediment concentration of 200
mg/L is unreasonable, this case was run to examine the effect of
sediment concentrations on radionuclide transport. Deep-ocean
sediment concentrations near the 3800-m site ranged from 0.02 to
0.1 mg/L, but there were also measurements indicating that
concentrations near the bottom could be as high as 100 mg/L.
Figures 5-34 to 5-36 show the computed radionuclide results.
Sediment-sorbed radionuclide concentrations per unit weight of
sediment are almost identical to those of Case 4. Thus, the
sediment sorbed concentration per unit volume of water in this case
is 1000 times higher than in Case 4. Although dissolved
concentrations calculated for Case 5 are smaller than those of Case
4, additional adsorption by the greater amount of suspended
sediment for Case 5 is still too small to be discernable in these
figures. These results indicate the predictable effect of
suspended sediment concentrations on the radionuclide distribution.
5.2.4 Simulation Results of Nonstratified Case 6
This is the same as Case 5 except that there is no
stratification caused by water temperature. Computed results are
shown in Figures 5-37 to 5-39, displaying the almost identical
distributions as those shown in Case 5. Close comparisons of
computed radionuclide concentrations of Cases 5 and 6 revealed up
to several percent differences in concentrations between these two
cases, due to slightly different eddy viscosity values calculated
for stratified and non-stratified cases. If we run a case with a
much smaller eddy viscosity calculated by stratified Case 3, we
expect to see a much more drastic difference in computed
radionuclides. Although under this study we did not apply the
FLESCOT code to three-dimensional cases, the model has been applied
to calculate these quantities in three-dimensional
shallowestuarine/coastal waters (Onishi et al. 1985).
5.10
-------
3000.00
o
c
o
.*-•
OT
Q
"5
o
0.000
PLOT FOR TYME =
173000.000
IxlCf10 m2/s2
o
3
83
o
u
0.000
Longitudinal Distance (m)
FIGURE 5.1. Calculated Turbulent Kinetic Energy Over the Whole Flow
Region for Nonstratified Case 1
B X
Ml
m
12000.000
-------
3000.00
V
O
c
O
n
b
"o
o
en
•
t—•
PO
0.000
PLOT FOR TYME =
173000.000
CO
CM
E
CM
t—I
O
i-H
X
tn
i-H
o
i—i
x
o
T)
O
•o
o
ID
O>
0.000
FIGURE 5.2.
Longitudinal Distance (m)
Calculated Turbulent Energy Dissipation Over the Whole
Flow Region for the Nonstratlfied Case 1
12000.000
-------
3000.00
0)
O
c
O
-*-»
CO
Q
"5
o
in
•
t-»
CO
0.000
PLOT FOR TYME =
173000.000
o
o
4)
IA
I
-------
68.460
PLOT FOR TYME =
173000.000
4)
U
c
o
-t->
OT
Q
"6
o
Ul
0.000
o
0.000
Longitudinal Distance (m)
TJ
O
•4^
O
o
!
12000.000
FIGURE 5.4. Calculated Turbulent Kinetic Energy Near the Bottom for
the Nonstratified Case 1
-------
68.460
PLOT FOR TYME =
173000.000
o>
u
c
o
•*-•
w
Q
8
01
•
>^
in
0.000
o
-Q
C
sx
0>
-------
68.460
PLOT FOR TYME
173000.000
4)
O
C
O
-f-»
(0
Q
"5
O
ITI
•
H-*
01
0.000
0.000
Longitudinal Distance (m)
S • -o
25?
12000.000
FIGURE 5.6. Calculated Eddy Viscosity Near the Bottom for the
Nonstratified Case 1
-------
en
Turbulent Kinetic Energy (m2/sec2) r2dla
o o run nr TDK -
o.ooo
r2dla
f
0.00 II.JS 22.30 33.75 45.00
VartLcaL DLstancs (m) r2dla
67.S3
FIGURE 5.7. Vertical Distribution of Calculated Turbulent Kinetic
Energy Near the Bottom for the Nonstratified Case 1
-------
Turbulent Energy Dlsslpltatlon (m2/sec
in
•
h-
00
O 8 FUR
I
1 I
t -
0.000 0.375
a.am
0.7SO 1.135 1.500
Vertical Distance (•)
1.87S 2.3
Figure 5.8.
Vertical Distribution of Calculated Turbulent Energy Dissipation Near the
Bottom for the Nonstratified Case 1
-------
Eddy viscosity r2dla
in
•
i-*
to
a run ta TIIC -
I
"8
8
0.00 11.25
0.000
SJ.7S 45.00
VartlcaL OLalonca (•)
SB.25 67.SO 78
FIGURE 5.9. Vertical Eddy Viscosity Near the Bottom for the
Nonstratified Case 1
-------
171
•
ro
o
68.460
PLOT FOR TYME
173000.000
o>
o
c
o
Q
"5
o
0.000
>>>>>>>>>>>>
> > » » >
0.000
Longitudinal Distance (m)
FIGURE 5.10. Calculated Velocity Distribution Near the Bottom for
the Nonstratified Case 1
o
=>
£
o
12000.000
-------
3000.00
O
O
-*••
tn
Q
"o
o
ro
0.000
PLOT FOR TYME =
173000.000
o
3
c <
o
o
0.000
Longitudinal Distance (m)
12000.000
FIGURE 5.11. Calculated Turbulent Kinetic Energy Over the Whole Flow
Region for Stratified Case 2
-------
3000.00
PLOT FOR TYME
173000.000
8
o
-*->
c
- o S
u
o
o
0.000
Longitudinal Distance (m)
12000.000
FIGURE 5.12. Calculated Turbulent Energy Dissipation Over the Whole
Flow Region for the Stratified Case 2
-------
3000.00
PLOT FOR TYME =
173000.000
O
O
-t->
W
Q
"o
O
ro
w
0.000
I/I
•o
O
O
41
VI
I
-------
68.460
PLOT FOR TYME =
173000.000
-------
68.46O
PLOT FOR TYME
173000.000
O
c
O
-•-•
CO
O
"5
O
ro
in
0.000
o
3
OO
TJ
O
O
0.000
Longitudinal Distance (m)
-------
68.460
PLOT FOR TYME
173000.000
O
I
•o
•
-t*
o
o
0.000
Longitudinal Distance (m)
12000.000
FIGURE 5.16. Calculated Eddy Viscosity Near the Bottom for the
Stratified Case 2
-------
ro
Turbulent Kinetic Energy (m2/sec2) r2d2
HTTUC -
0.000
r2d2
w fl
o "
I
So
o
-J
o.oo u.as aa.sn SI.TS «.oo
VertLcal DUtance (•) r2d2
SR.Z
67.50
FIGURE 5.17. Vertical Distribution of Calculated Turbulent Kinetic
Energy Near the Bottom for the Stratified Case 2
-------
in
•
PO
00
Turbulent Energy Dissipation (m2/sec3)
o s Pun nr TDK -
.
"8
~ B
§
* N
a. a
• 3
S N
0.000
r2d3
0.000 0.375 0.750 1.125 1.500
Vertical Distance (m) r2d3
1.87S
2.250
FIGURE 5.18.
Vertical Distribution of Calculated Turbulent Energy Dissipation Near the
Bottom for the Stratified Case 2
-------
in
Eddy vLscosLtij r2d2
o PUT m roc -
~ R
P
k.
0.00
11.
o.o
29.50 33.75 4S.OO
Vertical Distance (•)
sr.sa
n
FIGURE 5.19. Vertical Eddy Viscosity Near the Bottom for the
Stratified Case 2
-------
68.460
PLOT FOR TYME = 173000.000
o
c
o
•*»
in
Q
"o
o
en
g
0.000
0.000
Longitudinal Distance (m)
FIGURE 5.20. Calculated Velocity Distribution Near the
the Stratified Case 2
-------
3000.00
0)
O
O
-»-•
(0
b
"5
o
en
•
CO
0.000
PLOT FOR TYME =
173000.000
Ixli
0.000
FIGURE 5.21
Longitudinal Distance (m)
Calculated Turbulent Kinetic Energy Over the
Region for Stratified Case 3
-------
3000.00
o
c
o
-•-•
OT
Q
"6
o
en
•
CO
ro
0.000
PLOT FOR TYME =
173000.000
C\J
o
•-«
X
a
x
ro
r-l
I
O
t-H
X
o.ooo
Longitudinal Distance (m)
FIGURE 5.22.
Calculated Turbulent Energy Dissipation Over
Flow Region for the Stratified Case 3
-------
3000.00
PLOT FOR TYME =
173000.000
O
c
O
-•-»
-------
68.460
PLOT FOR TYME =
173000.000
o>
o
c
O
-»-•
(0
Q
"o
O
en
0.000
0.000
FIGURE 5.24.
Longitudinal Distance (m)
Calculated Turbulent Kinetic Energy Near th
the Stratified Case 3
-------
68.460
PLOT FOR TYME =
173000.000
I
o
2
Q
"5
o
en
0.000
0.000
FIGURE 5.25.
Longitudinal Distance (m)
Calculated Turbulent Energy Dissipation Ne<
for the Stratified Case 3
-------
68.460
PLOT FOR TYME =
173000.000
CO
en
O
C
O
-*-•
OT
o
O
o.ooo
o.ooo
Longitudinal Distance (m)
FIGURE 5.26. Calculated Eddy Viscosity Near the Bott<
Stratified Case 3
-------
en
•
to
Turbulent Kinetic Energy (m2/sec2) r
nr TOC -
r2d3
0.000 0.375 0.7SO I.IZ l.SOO 1.875
Vertical Distance (•) r2d3
FIGURE 5.27. Vertical Distribution of Calculated Turbul
Energy Near the Bottom for the Stratified
-------
run RT Tire -
o.ooo
I
1
S *J
1
1 I
vi
•
CO
00
1
0.000 0.375
0.750 1.125 1.500
Vertical Distance (m)
3.35
FIGURE 5.28.
Vertical Distribution of Calculated Turbulent Energy
the Bottom for the Stratified Case 3
-------
not w TOC '
O.DDO
tfl
<*>
to
JT
3 =
0.0
J.S
5.0 7.5 10.0
Vertical Distance (•)
12.5
FIGURE 5.29. Vertical Eddy Viscosity Near the Bottoi
Stratified Case 3
-------
68.480
PLOT FOR TYME =
173000.000
171
•
4k
O
O
a
-*->
(0
a
"o
o
0.000
0.000
FIGURE 5.30.
Longitudinal Distance (m)
Calculated Velocity Distribution Near th
the Stratified Case 3
-------
38.480
PLOT FOR TYME =
21600.000
10 pCi/L
8
o
"o
0.000
0.000
horizontal distance (m)
FIGURE 5.31.
Calculated Dissolved Radionurlide (
(in pCi/L) Near the Bottom for the
-------
38.460
PLOT FOR TYME
21600.000
IxlO"6 pCi/g
O1
•
•P»
ts>
o>
o
o
.*-»
CO
'•6
"o
o
0.000
0.000
horizontal distance (m)
FIGURE 5.32.
Calculated Suspended-SiIt-Sorbed Radionuc
(in pCi/g) Near the Bottom for the Strati
-------
38 460
PLOT FOR TYME =
216OO.OOO
VI
8
o
•*J
M
8
4>
0.000
0.000
5x!0"6 pCi/g
horizontal distance (m)
FIGURE 5.33. Calcualted Suspended-Clay-Sorbed Radionuc
(in pCi/g) Near the Bottom for the Strati
-------
38.460
PLOT FOR TYME =
21600.000
10 pCi/L
100
in
•
•p.
.52
TJ
"o
o
0.000
0.000
horizontal distance (m)
FIGURE 5.34.
Calculated Dissolved Radionuclide Con
(in pCi/L) Near the Bottom for the St
-------
38.460
PLOT FOR TYME
21600.0OO
Ixio"6 pCi/g
•pk
en
o>
o
o
-*-*
CO
o
o
0.000
0.000
FIGURE 5.35.
horizontal distance (m)
Calculated Suspended-SiIt-Sorbed Radionuc
(in pCi/g) Near the Bottom for the Strati
-------
38.460
PLOT FOR TYME
21600.000
5xlO~6 pCi/g
Ol
8
o
J2
-o
"o
4)
0.000
0.000
FIGURE 5.36.
horizontal distance (m)
Calculated Suspended-Clay-Sorbed Radionuclid
(in pCi/g) Near the Bottom for the Strati fie
-------
38.460
PLOT FOR TYME =
21600.000
10 pCi/L
u
o
-t->
in
T3
§
i.
4>
0.000
100
0.000
horizontal distance (m)
FIGURE 5.37. Calculated Dissolved Radionuclide (
(in pCi/L) Near the Bottom for the
-------
38.460
PLOT FOR TYME =
21600.000
lxlO~6 pCi/g
in
•
-P«
00
§
o
"o
O
V
0.000
0.000
FIGURE 5.38.
horizontal distance (m)
Calculated Suspended-SiIt-Sorbed Radionuclidt
(in pCi/g) Near the Bottom for the Stratifie<
-------
JS.46O
PLOT FOR TYME
216OO.OOO
4)
u
c
D
-t-»
CO
8
'•e
0)
0.000
bxio"6 pCi/g
0.000
horizontal distance (m)
FIGURE 5.39.
Calculated Suspended-Clay-Sorbed Radionucli
(in pCi/g) Near the Bottom for the Stratifi
-------
6. SUMMARY AND CONCLUSIONS
The TEMPEST/FLESCOT, FLOWER, Blumberg's, and RMA 10 models
would be appropriate candidates for regional modeling. Among these
models, TEMPEST/FLESCOT is currently the only model that solves
simultaneously for distributions of flow, turbulence (with the k-e
model), salinity, water temperature, sediment, dissolved
contaminants (including radionuclides, pesticides, heavy metals)
and sediment-sorbed contaminants.
Turbulence modeling has developed rapidly with the advent of
high-speed computers, but there are still limits on the complexity
of calculations that can be handled. Solving the Navier-Stokes
equations using higher order correlations requires a large number
of calculations and is not practical for modeling large areas where
numerous grid points must be calculated. Therefore, assumptions
must be made to solve the equations, which is where turbulence
modeling comes into play.
In turbulence modeling, most of the work has gone into trying
to define the turbulence correlations (Reynolds stresses) uf~uj and
U| 0. Once these terms are obtained, it is possible to solve the
equations for the flow field. The method used here for modeling
the turbulence correlations is the eddy viscosity/diffusivity
concept developed by Boussinesq in 1877. Zero, one, or two partial
differential equations must be solved to define the eddy
viscosity/diffusivity, depending on the assumptions made.
Zero-equation models do not require the solution of any
additional partial differential equation. The eddy
viscosity/diffusivity is defined by empirical means. One of the
most common methods being used is that of constant eddy
viscosity/diffusivity. This method is crude, especially in a deep
ocean model, because it is known that eddy viscosity/diffusivity
will vary across a flow field. A simple method for solving this
6.1
-------
problem is to specify values of eddy viscosity/diffusivity for
regions of the flow field, but this is still not very accurate.
Defining eddy viscosity/diffusivity as a function of shear stress
and depth is probably the simplest method that will give a
reasonably accurate representation of what is occurring in the flow
field. This method has been used in boundary layers, but there
will have to be some modification before it can accurately
represent what is occurring throughout the ocean.
One-equation models are the first to take into account the
history and transport effects in a flow field. This allows
turbulence to travel downstream, which was not possible when the
eddy viscosity/diffusivity was simply specified. The turbulent
kinetic energy (TKE) equation is used to define the velocity scale
in the eddy viscosity, and the eddy diffusivity is a function of
eddy viscosity and the Prandtl or Schmidt number. The main problem
with this method is that the length scale must still be defined
before a value for eddy viscosity can be obtained. In regions
where the length scale is clearly definable, as in the boundary
layer, this method works rather well. In regions where there is
no clearly defined length scale, empirical formulas have been
developed, but because rather extensive computational time is
required for these formulas, it is sometimes easier to move up to
the more complex two-equation models.
Two-equation models have partial differential equations for
both the velocity scale and the length scale. The k-e model, where
k is the same TKE equation developed to represent the one-equation
model and e is the rate of dissipation used to represent the length
scale, is probably one of the most common of this type. This
method will include historical and transport effects of turbulence.
Although the buoyancy term is built into the k-e model, there may
be potential problems in that it may not account fully for the
effects of stratification on the production/dissipation of the TKE.
6.2
-------
This problem, if it exists, may be solved by using an extended k-e
model with an additional buoyant-effect consideration that will
more accurately define the effects of stratification at the expens^
of computer time.
Because the two-equation models have the highest potential to
properly calculate the eddy viscosity/diffusion coefficients, we
have selected the most common two-equation model, the k-e model,
to calculate eddy viscosity in the deep ocean.
We have applied our model, TEMPEST/FLESCOT, with the k-e model
to a two-dimensional water body having a depth of 3000 m. The
preliminary model results demonstrate the adequacy of the k-e
model. Those results show that, near the bottom, vertical eddy
viscosity increases the most linearly with distance from the
bottom, and then decreases rapidly when the velocity distribution
becomes almost uniform vertically. Elsewhere, vertical eddy
viscosity calculated by the model is uniform and somewhat smaller
than the molecular viscosity due to upstream boundary conditions.
Comparisons among stratified and nonstratified cases show the
reducing effects of the buoyancy term on TKE and eddy viscosity,
resulting in a more reasonable level of calculated eddy viscosity
under stratified conditions in a deep-ocean environment.
The computer code also successfully predicted, at least
qualitatively, the leaching of radionuclides from a hypothetical
contamination source on the ocean bottom to a "clean" overlying
water-column. The code also predicted that some leached
radionuclides in the water-column were adsorbed by suspended silt
and clay, based on distribution coefficient and mass transfer
rates, while transported downstream. Furthermore, originally
"clean" bottom sediments, located downstream of contaminated bottom
sediments, were predicted to adsorb a fraction of the leached
radionuclides, thereby spreading the potential for long-term
sources of contamination.
6.3
-------
7. REFERENCES
Blumberg, A. F. 1977. "Numerical Model of Estuarine
Circulation." J. Hydraul. Div. ASCE HY3:295-311.
Blumberg, A. F, and H. J. Herring. 1986. Circulation
Modelling Using Curvilinear Coordinates. Report No.
81, Dyanalysis, Princeton,
New Jersey.
Bradshaw, P., ed. 1978. Turbulence. 2nd edition.
Bryan, K. 1969. "A Numerical Method for the Study
of the Circulation of the World Ocean." J. Coroput.
Phys. 4:347-376.
Davies, A. M. 1980. "Application of the Galerkin
Method to the Formulation of a Three-Dimensional
Nonlinear Hydrodynamic Sea Model." Appl. Math.
Modeling 4:425-455.
Ditoro, D. M., J. J. Fitzpatrick, and R. B. Thomann.
1981. Documentation for Water Quality Analysis
Simulation Program (WASP) and Model Verification
Program (MVP). Prepared for the Environmental
Research Laboratory, U.S. Environmental Protection
Agency, Duluth, Minnesota, by Hydroscience, Inc.,
Westwood, New Jersey.
Eraslan, A. H. , W. L. Lin, and R. D. Sharp. 1983.
FLOWER; A Computer Code for Simulating
Three-Dimensional Flow. Temperature, and Salinity
Conditions in Rivers. Estuaries. and Oceans.
NUREG/CR-3172, U.S. Nuclear Regulatory Commission,
Washington, D.C.
Freitas, C. J., A. N. Findikakis, and R. L. Street.
1985. "Three-Dimensional Simulation of a
Buoyancy-Driven Cavity Flow." In Proceedings
of the 1985 ASCE Hydraulics Division Specialty
Conference. Lake Buena Vista, Florida, August 12-17,
1985, pp. 1101-1106.
Gawain, T. H. , and J. W. Pritchett. 1970. "A Unified
Heuristic Model of Fluid Turbulence." J. Comput. Phvs.
5:383-405.
7.1
-------
Geyer, W. R. 1973. The Time-Dependent Dynamics of
a Salt Wedge. Special Report 101, College of
Oceanography and Fishery Sciences. University of
Washington, Seattle, Washington.
Grant, W. D., and O. S. Madsen. 1979. "Combined
Waves and Current Interaction with a Rough Bottom."
J. of Geophvs. Res. 84(4)1797-1808.
Holland, W. R., and L. B. Liu. 1975. "On the
Generation of Mesoscale Eddies and Their Contribution
to the Oceanic General Circulation." J. Phys.
Oceanogr. 5(4):642-669.
Jones, W. P., and B. E. Launder. 1973. "The
Calculation of Low-Reynolds-Number Phenomena with a
Two-Equation Model of Turbulence." Int. J. Heat Mass
Transfer 16:1119-1130.
King, I. P. 1982. "A Three Dimensional Finite
Element Model for Stratified Flow." In Finite Element
Flow Analysis. T. Kawai, ed. , pp. 513-527. University
of Tokyo Press, Japan.
Koplik, C. M., M. F. Kaplan, J. Y. Nalbandian, J. H.
Simonson, and P. G. Clark. 1984. User's Guide to
MARINRAD; Model for Assessing the Consequences of
Release of Radioactive Material into the Oceans.
SAND83-7104. Prepared for Sandia National
Laboratories, Albuquerque, New Mexico, by Analytic
Sciences Corp., Reading, Massachusetts.
Krone, R. B. 1962. Flume Studies of the Transport
of Sediment in Estuarial Shoaling Processes.
Hydraulic Engineering Laboratory and Sanitary
Engineering Research Laboratory, University of
California at Berkeley, Berkeley, California.
Leenderste, J. J. , and S. K. Liu. 1975. "A
Three-Dimensional Model for Estuaries and Coastal
Seas." In Aspects of Computation,. Vol. n.
R-1784-OWRT, RAND Corp., Santa Monica, California.
Mellor, G. L., and T. Yamada. 1974. "A Hierarchy of
Turbulence Closure Models for Planetary Boundary
Layers." J. Atmos. Sci. 31:1791-1806.
7.2
-------
Mellor, G. L., and T. Yamada. 1982. "Development of
a Turbulence Closure Model for Geophysical Fluid
Problems." Rev. Geophys. Space Phys. 20(4):851-875.
Munk, W. J., and E. R. Anderson. 1948. "Note on the
Theory of the Thermocline." J. Mar. Res.
7(3):276-295.
Norton, W. R. , and I. P- King. 1977. Operating
Instructions for the Computer Program RMA2. RMA 6290,
Resource Management Associates,
Lafayette, California.
Onishi, Y. 1981. "Sediment and Contaminant Transport
Model." J. Hvdraul. Div. ASCE 107 (HY9):1089-1107.
Onishi, Y., and D. S. Trent. 1982. Mathematical
Simulation of Sediment and Radionuclide Transport in
Estuaries. NUREG/CR-2423. Prepared by Pacific
Northwest Laboratory for U.S. Nuclear Regulatory
Commission, Washington, D.C.
Onishi, Y., D. S. Trent, and A. S. Koontz. 1985a.
"Three-Dimensional Simulation of Flow and Sewage
Effluent Migration in the Strait of Juan De Fuca,
Washington." In Proceedings of the 1985 ASCE
Environmental Engineering Conference. Boston,
Massachusetts, pp. 1006-1013.
Onishi, Y., R. M. Ecker, D. S. Trent, S. B. Yabusaki
and L. W. Vail. 1985b. Oceanographic Modeling of the
Beaufort Sea in Relation to the Proposed Lisburne
Development. Prepared for ARCO Alaska, Inc., by
Battelle, Pacific Northwest Laboratories, Richland,
Washington.
Onishi, Y., and D. S. Trent. 1985. Three-Dimensional
Simulation of Flow. Salinity. Sediment. and
Radionuclide Movements in the Hudson River Estuary.
Proceedings of the ASCE Hydraulics Division Conference
at Lake Buena Vista, Florida, August 1985, pp. 1095-
1100.
Onishi, Y., L. F. Hibler, and C. R. Sherwood. 1987.
Review of Hydrodynamic and Transport Models and Data
Collected near the Mid-Atlantic Low-Level Radioactive
Waste Disposal Sites. PNL-6331, Pacific Northwest
Laboratory, Richland, Washington.
7.3
-------
Partheniades, E. 1962. A study of Erosion and
Deposition of Cohesive Soils in Salt Water.
Dissertation, University of California, Berkeley,
California.
Robinson, A. R., and M. G. Marietta. 1985. Research,
Progress and the Mark A Box Model for Physical.
Biological and Chemical Transports. SAND84-0646,
Sandia National Laboratories, Albuquerque, New Mexico.
202 p.
Rodi, W. 1984. Turbulence Models and Their
Application in Hydraulics. University of Karlsruhe,
Karlsruhe, Federated German Republic.
Rodi, W. 1987. "Examples of Calculation Methods for
Flow and Mixing in Stratified Fluids." J. Geophvs.
Res. 92(C5):5305-5328.
Sheng, Y. P., W. Lick, R. T. Gedney, and F. B. Molls.
1978. "Numerical Computation of a Three-Dimensional
Circulation in Lake Erie: A Comparison of a
Free-Surface Model and a Rigid-Lid Model." J. Phys.
Oceanogr. 8:713-727.
Shepard, J. G. 1978. "A Simple Model for the
Dispersion of Radioactive Wastes Dumped on the
Deep-Sea Bed." Marine Sci. Comm. 4(4):293-327.
Simons, T. J. 1973. "Development of
Three-Dimensional Numerical Models of the Great
Lakes." Scientific Series No. 12, Canada Center for
Inland Waters, Burlington, Ontario, Canada.
Smith, J. H., W. R. Mabey, N. Bohonos, B. R. Holt, S.
S. Lee, T. W. Chou, D. C. Bomberger, and T. Mill.
1977. Environmental Pathways of Selected Chemicals
in Freshwater Systems. Part l! Background and
Experimental Procedures. EPA-6000/7-77-H3. Prepared
by SRI International for the U.S. Environmental
Protection Agency, Washington, D.C.
Spaulding, M. L., and D. Parish. 1984. "A Three-
Dimensional Numerical Model of Particulate Transport
for Coastal Water." Continental Shelf T?OC 3(1)755-
67.
7.4
-------
Tee, K. T. 1981. "A Three-Dimensional Model for
Tidal and Residual Currents in Bays." In Transport
Models for Inland and Coastal Waters. H. B. Fischer,
ed., pp. 284-309, Academic Press, New York.
Trent, D. S., L. L. Eyler, and M. J- Budden. 1989.
TEMPEST - A Three-Dimensional Time-Dependent Computer
Program for Hydrothermal Analysis. PNL-4348, Vol. 1,
Rev- 2, Pacific Northwest Laboratory, Richland,
Washington.
U.S. Atomic Energy Commission. 1975. Overall Safety
Manual. Vol. 2. Space Nuclear Systems Division,
Washington, D.C.
Vanoni, V- A. 1975. Sedimentation Engineering.
Prepared by the ASCE Task Committee for the
Preparation of the Manual on Sedimentation,
Sedimentation Committee of the Hydraulics Division,
ASCE, New York.
Walker, H. A., J. A. Nacito, J. F. Paul,
V. J. Bierman, and J. H. Gentile. 1985. Methods for
Waste Load Allocation of Municipal Sewage Sludge at
the 106-Mile Ocean Disposal Site. Environmental
Research Laboratory, U.S. Environmental Protection
Agency, Narragansett, Rhode Island.
Wang, J. D., and J. J. Connor. 1975. Mathematical
Modeling of Near Coastal Circulation. Report No. 200,
Department of Civil Engineering, Massachusetts
Institute of Technology, Cambridge, Massachusetts.
Wang, D. P. 1982. "Development of a
Three-Dimensional, Limited-Area (Island) Shelf
Circulation Model." J. Phys. Oceanogr. 12:605-617.
Weatherly, G. L., and P. J. Martin. 1978. "On the
Structure and Dynamics of the Oceanic Bottom Boundary
Layer." J. Phys. Oceanogr. 8:557-570.
Webb, G. A. M. , and F. Morley. 1973. A Model for the
Evaluation of Deep Ocean Disposal of Radioactive
Waste. National Radiological Protection Board, p. 23.
7.5
-------