UFL/COEL92/020
PHYSICAL EFFECTS OF VEGETATION ON WINDBLOWN SAND
IN THE COASTAL ENVIRONMENTS OF FLORIDA
By
D. M. Sheppard
Department of Coastal and Oceanographic Engineering
and
A. W. Niedoroda
Environmental Science and Engineering Inc.
February, 1992
Physical Effects of Vegetation on WindBlown Sand
in the Coastal Environments of Florida
D.M. Sheppard
Department of Coastal and Oceanographic Engineering
University of Florida
A.W. Niedoroda
Environmental Science and Engineering Inc.
February, 1992
Form Approved
REPORT DOCUMENTATION PAGE OM No. 70o0188
PuoKC reortiNq Burden for utr colw ton of information etimateo to aver e I houw or reasons. including Xe tmte for reveewing imtrucon searching extstlnq data sourfe.
gathnermq and mafitaining the data needed. and compltia and revewjnn the collection Of information. Send comments reIgrdtnr thrl s de estimate or any other as o f Ihn
cotlection of information. indudhng suggestint for edduing tl burden. to Wa'shnqtorl MHedlourterm Sfrcte. O*sctorate for Imformatron Oeratomr and Relort 1215I Jefferson
DOav Hofgtn. Suete 1204. Arflgton. VA 222024302. and to the OfKte of Mana ent and Budget. P0aperot Redwctin o l eot (070440188). Wahington. OC 200)3.
1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 13. REPORT TYPE AND DATES COVERED
February 1992 Final Report
4. TITLE AND SUBTITLE 5. FUNDING NUMBERS
Physical Effects of Vegetation on WindBlown Sand C No. C6851
in the Coastal Environments of Florida
6. AUTHORS)
D.M. Sheppard
A.W. Niedoroda
7. PERFORMING ORGANIZATION NAME(S) AND AODRESS(ES) 8. PERFORMING ORGANIZATION
REPORT NUMBER
*University of Florida
Gainesville, Florida 32611 PR 4910451123312
*Environmental Science & Engineering Inc. UPN No. 90100411
Gainesville, Florida
9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESSES) 10. SPONSORING /MONITORING
AGENCY REPORT NUMBER
U.S. Dept. of Commerce/NOAA Dept. of Eng. Reg. AGEN REPORTNUMBER
OCRM Coastal Management
1825 Connecticut A., N.W. 2600 Blair Stone Rd.
Washington, D.C. 20235 Tallahassee, FL 32399
11. SUPPLEMENTARY NOTES
12a. DSTRIBUTIONIAVAILA8IUTY STATEMENT 12b. DISTRIBUTION CODE
13. ABSTRACT (Maximum 200wods)
One and twodimensional, second order turbulence plant canopy flow
models were developed for the purpose of estimating the effect of
coastal vegetation on wind blown sand transport. The computer program
that solves the governing differential equations uses measured leaf
area density profiles and drag coefficients for crop plants similar in
shape and size to the more common coastal vegetation in Florida.
14. SUBJECT TERMS 15. NUMBER OF PAGES
Eolian sand transport / wind blown sand transport / 54
16. PRICE CODE
canopy flow / sand trapping.
17. SECURITY CLASSIFICATION 18. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 20. LIMITATION OF ABSTRACT
OF REPORT OF THIS PAGE OF ABSTRACT
Standard Form 298 (Rev. 289)
Prlcrsbed by ANSI SId. Z3918
298_102
SN N 7540012805 0
I
Partial funds for this project were provided by the Department of Environmental Reg
ulation, Office of Coastal Management using funds made available through the National
Oceanographic and Atmospheric Administration under the Coastal Zone Management Act
of 1972, as amended.
1
Final Report
Physical Effects of Vegetation on WindBlown Sand
in the Coastal Environments of Florida
Principal Investigator: D.M. Sheppard
Department of Coastal and Oceanographic Engineering
University of Florida
CoPrincipal Investigator: A.W. Niedoroda
Environmental Science and Engineering Inc.
UF Project No. 4910451123312
UPN No. 90100411
Contract No. C 6851
Department of Coastal and Oceanographic Engineering
University of Florida
Gainesville, Florida 32611
February, 1992
Executive Summary
This report presents the results of a study to quantify the effects of coastal dune vege
tation on trapping and retaining wind blown sand. One and two dimensional plant canopy
flow and sand transport models were developed. Both flow models include a second order
closure turbulence scheme based on the model developed by Wilson and Shaw (1977). The
onedimensional model allows a vertical distribution of leaf area density and a variable plant
spacing but assumes a fully developed wind flow field in the canopy. The twodimensional
model can account for variations in plant type and spacing in the direction of wind flow
and solves for the flow and sand transport in the transition regions at the beginning and
end of the canopy as well as the interior (fully developed) region. The onedimensional
model forms the core of the completely redesigned Florida Department of Natural Resources
(FDNR) computer porgram, VEG, developed under a previous contract. The new VEG pro
gram (Version 2.0) uses parameterized results from numerical experiments performed with
the twodimensional model thus some of the more important twodimensional features are
included. No experimental data for leaf area density profiles, drag coefficients, etc. exists
in the literature for dune vegetation. Data for crop plants similar in shape and size to the
more common Florida dune vegetation has been used in VEG version 2.0.
PHYSICAL EFFECTS OF VEGETATION ON
WINDBLOWN SAND IN THE COASTAL
ENVIRONMENTS OF FLORIDA
Introduction
The objective of this research is to provide a means for predicting the sand trapping and
retention characteristics of the more common species of dune vegetation. This stored sand
is an important natural resource which aids in preventing long term coastal erosion and
acts as a barrier to storm surge during extreme storm events. The initial phases of this
year's work began in November of 1990 under a previous DNR contract and continued
under the present contract. The products of the first study include: 1) establishing the
physical characteristics of vegetation that control nearground air flow and eolian sand
transport, 2) evaluating the shape characteristics of coastal dunes along the entire coast of
Florida, 3) establishing the dominant types of natural vegetation in Florida's coastal areas,
4) developing a methodology for computing the wind blown sand trapping and retention
properties of various types of dune vegetation, 5) estimating the physical characteristics
and aerodynamic parameters associated with each of the dominant natural plant types, 6)
establishing an approximate index of the root binding effectiveness of the plant types, 7)
writing a computer program for estimating the rate of sand accretion or erosion and the
root binding effectiveness for different heights and densities of natural coastal vegetation,
and 8) starting preliminary development work on a onedimensional numerical model of air
flow through plant canopies. The approach taken in the first year's work was new and
innovative and provided, for the first time, a quantitative way to address this problem.
The initial analysis of the canopy flow problem was analytical, i.e. it involved solving the
governing differential equations analytically. While this approach has certain advantages, it
has severe limitations when more complex vegetation, assemblages of vegetation and
ground relief must be considered. Numerical techniques must be used for these more
complex problems. This year's work addressed some of these limitations. The terms model,
algorithm and (computer) program are often used interchangeably. In this report the term
model will be used to designate a "mathematical model" of a process, algorithm will be
used in connection with a (usually numerical) scheme for solving an equation and
(computer) program will denote the computer code used to implement the numerical
algorithms, solve the governing equations and present the output.
Objectives and Solution Methods
The main objectives of the present study were to expand the usefulness of the 1990 model
and computer program by 1) allowing more realistic plant leaf area distributions and plant
assemblages and spacings and 2) accounting for flow transitional effects by extending the
model and computer program to two spatial dimensions. The first objective was achieved
by completing the development of the onedimensional canopy flow and eolian sand
transport model and program. Addition of the second order turbulence closure model
greatly improved the accuracy of the computed flow field and sand transport and allowed
more realistic representations of the plants, combinations of plants and plant spacings. To
satisfy the second objective, a twodimensional model and program was developed using an
extension of the turbulence algorithm used in the onedimensional case. The two
dimensional analysis is necessary since the canopy depth (length in the direction of the
wind) to plant height ratio is small for many cases of interest to FDNR engineers. Under
these conditions the flow and sand transport are always in transition and never reach the
equilibrium state required for use of the onedimensional model. A brief description of the
capabilities and limitations of the one and twodimensional models and programs is given
below. The governing differential equations and the algorithms used in their solution are
given in the appendices. A case study using the programs developed in this study is
presented in Appendix E.
Model Capabilities and Limitations
The onedimensional model and computer program serves two purposes. First it is a stand
alone tool for the computation of wind blown sand accretion or erosion rates in dune
vegetation canopies whose depth is large compared to the plant height. Second it is a
desirable (if not necessary) step in the development of a twodimensional model. The
turbulence models used in both the one and twodimensional models are based on the work
of Wilson and Shaw (1977). The numerical solution of the governing equations for these
models allow variations in the distribution of leaf areas with height above the ground, allow
an assemblage of up to six different plants to be considered, and more accurately account
for plant spacing than the analytical solution.
For the onedimensional model the flow field is assumed to be fully developed, i.e.
independent of the coordinate in the direction of the wind. This is not a severe limitation
in the case of grass or small shrubs but can yield the model invalid for limited areas of tall
plants such as trees. The initial intent was to use the second order turbulence program to
perform numerical experiments. The results of these would then be parameterized and
incorporated into the existing VEG PROGRAM. During the course of the study it became
obvious that rewriting the VEG PROGRAM to include the onedimensional second order
Objectives and Solution Methods
The main objectives of the present study were to expand the usefulness of the 1990 model
and computer program by 1) allowing more realistic plant leaf area distributions and plant
assemblages and spacings and 2) accounting for flow transitional effects by extending the
model and computer program to two spatial dimensions. The first objective was achieved
by completing the development of the onedimensional canopy flow and eolian sand
transport model and program. Addition of the second order turbulence closure model
greatly improved the accuracy of the computed flow field and sand transport and allowed
more realistic representations of the plants, combinations of plants and plant spacings. To
satisfy the second objective, a twodimensional model and program was developed using an
extension of the turbulence algorithm used in the onedimensional case. The two
dimensional analysis is necessary since the canopy depth (length in the direction of the
wind) to plant height ratio is small for many cases of interest to FDNR engineers. Under
these conditions the flow and sand transport are always in transition and never reach the
equilibrium state required for use of the onedimensional model. A brief description of the
capabilities and limitations of the one and twodimensional models and programs is given
below. The governing differential equations and the algorithms used in their solution are
given in the appendices. A case study using the programs developed in this study is
presented in Appendix E.
Model Capabilities and Limitations
The onedimensional model and computer program serves two purposes. First it is a stand
alone tool for the computation of wind blown sand accretion or erosion rates in dune
vegetation canopies whose depth is large compared to the plant height. Second it is a
desirable (if not necessary) step in the development of a twodimensional model. The
turbulence models used in both the one and twodimensional models are based on the work
of Wilson and Shaw (1977). The numerical solution of the governing equations for these
models allow variations in the distribution of leaf areas with height above the ground, allow
an assemblage of up to six different plants to be considered, and more accurately account
for plant spacing than the analytical solution.
For the onedimensional model the flow field is assumed to be fully developed, i.e.
independent of the coordinate in the direction of the wind. This is not a severe limitation
in the case of grass or small shrubs but can yield the model invalid for limited areas of tall
plants such as trees. The initial intent was to use the second order turbulence program to
perform numerical experiments. The results of these would then be parameterized and
incorporated into the existing VEG PROGRAM. During the course of the study it became
obvious that rewriting the VEG PROGRAM to include the onedimensional second order
turbulence code would be highly advantageous. There are, however, significant differences
in a program to be used for numerical experiments by its developer and one that is used by
many users for a wide variety of input conditions. Ways to improve the speed and
dependability of the code were undertaken (i.e. steps were taken to be sure that the time
increments and grid size were such as to insure convergence for the range of input
parameters anticipated, etc.). A new onedimensional VEG PROGRAM called VEG
Version 1.5 was developed. Efforts were made to make the changes transparent to the user
of Version 1.0 and to maintain and enhance the ease at which the program can be used. A
number of the help screens now include graphics to assist with the input parameters and to
help interpret the output parameters.
The twodimensional program not only solves for the flow and sand transport in the
transitional zones but also provides the framework for inclusion of ground relief in future
models. Even though this model has been written with speed and efficiency in mind it still
requires approximately three hours of CPU time on the large University of Florida IBM
mainframe computer (IBM model 3090600J). Details of the twodimensional model,
governing equations and the numerical algorithms used in the solution are given in
Appendix C. The twodimensional program was used to perform numerical experiments
and the results parameterized for inclusion in an expanded version of the VEG
PROGRAM. This expanded version is called VEG Version 2.0. Some flexibility is lost
using this approach, but until (micro) computers become significantly faster it will not be
feasible to do otherwise. Parameterized results from the twodimensional solution include:
1) the length within the canopy over which sand is accumulated, 2) the rate of sand
accumulation in the canopy, 3) the distance (length in the direction of wind flow measured
from the front edge of the canopy) affected by the canopy, and 4) the length over which
sand is eroded and the rate of erosion down stream of the canopy (this assumes that the
surface down stream of the canopy is bare sand and free to erode). Knowledge of these
quantities can be extremely important to the engineer assessing the impact of a proposed
modification to the vegetation at a particular coastal site.
Work on the inclusion of an irregular or wavy boundary in the twodimensional model and
program was also initiated as part of this study. This is by no means a trivial task due to
the complexity of the solution scheme needed to solve for the pressure field for this case.
The first step, i.e. the solution to the laminar flow equations, was accomplished in this
study and some of the results are given in Appendix D. Since many beach dune systems
have sufficient ground relief that they cannot be approximated as flat, this work needs to
be taken to completion.
Finally, a case study of a hypothetical permit application using the programs developed in
this study is presented in Appendix E. It is somewhat idealized but illustrates some of the
ways these programs can be used as tools in the evaluation of vegetation modification
impact on the fragile beachdune system.
Summary and Recommendations for Future Work
The use of sophisticated computer models has become a necessary part of engineering
practice. As the speed and capacity of computers increase and as numerical algorithms
improve, the complexity of mathematical models that can be solved increases. However,
regardless of the complexity of the model it is still an idealization of the real system and
thus must 1) be verified with good experimental data and 2) be used with good engineering
judgement in the solution of problems. The models developed in this year's work greatly
enhance the capability of FDNR engineers' ability to estimate the sand trapping and
retention capabilities of various types of dune vegetation. Advantages of the new
onedimensional VEG PROGRAM (Version 1.5) over the earlier program (Version 1.0)
include: being able to analyze more realistic plant canopy distributions (LAD profiles);
being able to accurately analyze plant assemblages (up to six plants); being able to more
accurately account for plant spacing and; making it much easier to incorporate the
parameterized results from the twodimensional model. The bottom line difference is more
accurate results for situations where fully developed flow exists throughout most of the
canopy.
The second program, VEG Version 2.0, contains everything in Version 1.5 plus
parameterized results from numerical experiments with the twodimensional model.
Version 2.0 has not been tested to the extent that Version 1.5 has and should be considered
as experimental at this point. As stated above the twodimensional effects are particularly
important in situations where the vegetation height is large compared to the horizontal
depth of vegetation in the wind flow direction. Once these programs are thoroughly tested
with measured data, they can be used by FDNR engineers in their difficult and complex
task of evaluating requests for modifying and developing Florida's coastal areas.
Good laboratory and field data for some of the parameters used in, and predicted by, these
programs, such as plant drag coefficients, leaf area density distributions, transitional wind
velocity and turbulence intensities, are scarce or nonexistent. Measurement of these
quantities should be the next step in this research effort followed by or parallel with the
inclusion of ground relief in the twodimensional model.
Acknowledgements
Several people contributed to the success of this project. These include:
Dr. Christopher Reed, Post Doctoral Fellow and computational fluid dynamicist,
made significant contributions in the development of the mathematical models,
numerical algorithms, and computer code.
Summary and Recommendations for Future Work
The use of sophisticated computer models has become a necessary part of engineering
practice. As the speed and capacity of computers increase and as numerical algorithms
improve, the complexity of mathematical models that can be solved increases. However,
regardless of the complexity of the model it is still an idealization of the real system and
thus must 1) be verified with good experimental data and 2) be used with good engineering
judgement in the solution of problems. The models developed in this year's work greatly
enhance the capability of FDNR engineers' ability to estimate the sand trapping and
retention capabilities of various types of dune vegetation. Advantages of the new
onedimensional VEG PROGRAM (Version 1.5) over the earlier program (Version 1.0)
include: being able to analyze more realistic plant canopy distributions (LAD profiles);
being able to accurately analyze plant assemblages (up to six plants); being able to more
accurately account for plant spacing and; making it much easier to incorporate the
parameterized results from the twodimensional model. The bottom line difference is more
accurate results for situations where fully developed flow exists throughout most of the
canopy.
The second program, VEG Version 2.0, contains everything in Version 1.5 plus
parameterized results from numerical experiments with the twodimensional model.
Version 2.0 has not been tested to the extent that Version 1.5 has and should be considered
as experimental at this point. As stated above the twodimensional effects are particularly
important in situations where the vegetation height is large compared to the horizontal
depth of vegetation in the wind flow direction. Once these programs are thoroughly tested
with measured data, they can be used by FDNR engineers in their difficult and complex
task of evaluating requests for modifying and developing Florida's coastal areas.
Good laboratory and field data for some of the parameters used in, and predicted by, these
programs, such as plant drag coefficients, leaf area density distributions, transitional wind
velocity and turbulence intensities, are scarce or nonexistent. Measurement of these
quantities should be the next step in this research effort followed by or parallel with the
inclusion of ground relief in the twodimensional model.
Acknowledgements
Several people contributed to the success of this project. These include:
Dr. Christopher Reed, Post Doctoral Fellow and computational fluid dynamicist,
made significant contributions in the development of the mathematical models,
numerical algorithms, and computer code.
Ms. Tracy Burger, Master of Science student and Graduate Assistant.
Mr. Subarna Malakar, Associate Research Engineer and Head of the Computations
Laboratory in the Coastal and Oceanographic Engineering Department.
Ms. Erin F. Parker, Computer Systems Manager, Environmental Science
Engineering, Inc., Tampa office.
Ms. Laura Dickinson, secretary and word processing specialist, assisted with the
organization of the project and typed all reports and documents associated with the
project.
The authors are most grateful to these people for their hard work and dedication to the
project. Dr. Alfred Devereaux and his capable engineering staff provided valuable input
throughout the project. In addition to his technical input, Mr. Paden Woodruff also
contributed by serving as project monitor and administrator.
Appendix A
Dune Vegetation Bibliography
Appendix A
Dune Vegetation Bibliography
Anderson, D.E., S.B. Verma, R.J. Clement, D.D. Baldocchi and D.R. Matt, (1986)
"Turbulence Spectra of Co2, Water Vapor, Temperature and Velocity over a
Deciduous Forest." Agricultural and Forest Meteorology, V.38, 8199.
Anderson, D.A., J.C. Tannehill, R.H. Pletcher, (1984) Computational Fluid Mechanics
and Heat Transfer, Hemisphere Publishing Corporation, New York,
Avissar, R., M.D. Moran, G. Wu, R.N. Meroney and R.A. Pielke, (1990) "Operating
Ranges of Mesoscale Numerical Models and Meteorological Wind Tunnels for the
Simulation of Sea and Land Breezes." BoundaryLayer Meteorology, V.50, 227275.
Bagnold, R.A., (1941) Libyan Sands, Hodder and Stoughten, St. Paul's House, London,
E. C., V.4, 288 p.
Baldocchi, D.D. and B.A. Hutchison, (1988) "Turbulence in an Almond Orchard: Spatial
Variations in Spectra and Coherence." BoundaryLayer Meteorology, V.42, 293311.
Baldocchi, D.D. and T.P. Meyers, (1988) "A Spectral and LagCorrelation Analysis of
Turbulence in a Deciduous Forest Canopy." BoundaryLayer Meteorology, V.45,
3158.
Baldocchi, D.D. and T.P. Meyers, (1987) "Turbulence Structure in a Deciduous Forest."
BoundaryLayer Meteorology, V.43, 345364.
Baldocchi, D.D. and B.A. Hutchison, (1987) "Turbulence in an Almond Orchard: Vertical
Variations in Turbulent Statistics." BoundaryLayer Meteorology, V.40, 127146.
Baldocchi, D.D., S.B. Verma and N.J. Rosenberg, (1983) "Characteristics of Air Flow
Above and Within Soybean Canopies." BoundaryLayer Meteorology, V.25, 4354.
Baldocchi, D.D., (1982) "Mass and Energy Exchanges of Soybeans: MicroclimatePlant
Architectural Interactions." Dissert. U. of Nebraska, Lincoln.,
Briggs, W.L., (1987) A Multigrid Tutorial, Lancaster Press, Lancaster, Pennsylvania, 90p.
Busnaina, A.A., G. Ahmadi and S.J. Chowdhury, (1987) "Two Equation Turbulence
Model Consistent with the Second Law." AIAA Journal, V.25, No.12, 15431544.
Eghlima, A. and C. Kleinstreuer, (1985) "Numerical Analysis of Attached Turbulent
Boundary Layers along Strongly Curved Surfaces." AIAA Journal, V.23, No.2,
177184.
Fast, J.D. and E.S. Takle, (1988) "Evaluation of an Alternative Method for Numerically
Modeling Nonhydrostatic Flows over Irregular Terrain." BoundaryLayer Mete
orology, V.44, 181 206.
Fast, J.D. and E.S. Takle, (1988) "Application of A Quasi Nonhydrostatic
Parameterization for Numerically Modeling Neutral Flow over an Isolated Hill."
Boundary Layer Meteorology, V.44, 285304.
Finnigan, J.J. and P.J. Mulhearn, (1978) "Modelling Waving Crops in a Wind Tunnel."
BoundaryLayer Meteorology, V.14, 253 277.
Finnigan, J.J., (1979) "Turbulence in Waving Wheat." Boundary Layer Meteorology,
V.16, 181211.
Fortin, A., M. Fortin and J.J. Gervais, (1987) "A Numerical Simulation of the Transition
to Turbulence in a Two Dimensional Flow." J. of Computational Physics, V.70,
295310.
Garratt, J.R., (1990) "The Internal Boundary Layer A Review." BoundaryLayer Mete
orology, V.50, 171203.
Gong, W. and A. Ibbetson, (1989) "A Wind Tunnel Study of Turbulent Flow over Model
Hills." BoundaryLayer Meteorology, V.49, 113148.
Halldin, S., (1985) "Leaf and Bark Area Distribution in a Pine Forest." The
ForestAtmosphere Interaction, 3958.
Haworth, D.C. and S.B. Pope, (1986) "A Generalized Langevin Model for Turbulent
Flows." Phys. Fluids, V.29, No.2, 387405.
Hogstrom, U., (1988) "NonDimensional Wind and Temperature Profiles in the
Atmospheric Surface Layer: A ReEvaluation." BoundaryLayer Meteorology, V.42,
5578.
Horikawa, K., S. Hotta and S. Kubota, (1986) "Field Measurement of Vertical
Distribution of Wind Speed with Moving Sand on a Beach." Coastal Engineering in
Japan, V.29, 163178.
Horikawa, K., S. Hotta, S. Kubota and S. Katori, (1984) "Field Measurement of Brown
Sand Transport Rate by Trench Trap." Coastal Engineering in Japan, V.27, 213232.
Horikawa, K., S. Hotta, S. Kubota and S. Katori, (1983) "On the Sand Transport Rate by
Wind on a Beach." Coastal Engineering in Japan, V.26, 101120.
Horikawa, K., S. Hotta and S. Kubota, (1982) "Experimental Study of Blown Sand on a
Wetted Sand Surface." Coastal Engineering in Japan, V.25, 177195.
Hotta, S., N.C. Kraus and K. Horikawa, (1987) "Function of Sand Fences in Controlling
WindBlown Sand." Proc., Coastal Sediments '87, ASCE,, V.1, 1214.
Hotta, S., S. Kubota, S. Katori and K. Horikawa, (1984) "Sand Transport by Wind on a
Wet Sand Surface." Proc., Nineteenth Coastal Engineering Conf., V.2, 12651281.
Hutchison, B.A., et. al., (1986) "The Architecture of a Deciduous Forest Canopy in
Eastern Tennessee, U.S.A." J. of Ecology, V.74, 635646.
Inoue, K. and Z. Uchijima, (1979) "Experimental Study of Microstructure of Wind
Turbulence in Rice and Maize Canopies." Bull. Natl. Inst. Agric. Sci., Ser. A,
No.26, 188.
Isobe, S., (1972) "A Spectral Analysis of Turbulence in a Corn Canopy." Bull. Nat. Inst.
Agr. Sci., Ser.A, No.19., 101 112.
Kalma, J.D. and G. Stanhill, (1971) "The Climate of an Orange Orchard: Physical
Characteristics and Microclimate Relationships." Agricultural Meteorology, V.10,
185201.
Kalma, J.D., (1970) "Some Aspects of the Water Balance of an Irrigated Orange
Plantation." The Volcani Institute of Agricultural Research,
Kim, J. and P. Moins, (1985) "Application of a Fractional Step Method to Incompressible
NavierStokes Equations." J. of Computational Physics, V.59, 308323.
Kimura, F. and R. Manins, (1988) "Blocking in Periodic Valleys." Boundary Layer Mete
orology, V.44, 137169.
Kleinstreuer, C., A. Eghlima and J.E. Flaherty, (1982) "New HigherOrder
BoundaryLayer Equations for Laminar and Turbulent Flow Past Axisymmetric
Bodies." AFOSRTR820588, 33 p.
Landsberg, J.J. and G.B. James, (1971) "Wind Profiles in Plant Canopies: Studies on an
Analytical Model." J. of Applied Ecology, V.8, 729741.
Legg, B.J. and I.F. Long, (1975) "Turbulent Diffusion within a Wheat Canopy: II.
Results and Intrepretation." Quarterly J. of the Royal Meteorological Society, V.101,
No.429, 611628.
Li, Z., J.D. Lin and D.R. Miller, (1990) "Air Flow over and Through a Forest Edge: A
SteadyState Numerical Simulation." BoundaryLayer Meteorology, V.51, 179197.
Maitani, T., (1989) "WaveLike Wind Fluctuations Observed in the Stable Surface Layer
over a Plant Canopy." BoundaryLayer Meteorology, V.48, 1931.
Mellor, G.L., (1973) "Analytic Prediction of the Properties of Stratified Planetary Surface
Layers." J. of the Atmospheric Sciences, V.30, 10611069.
Moeng, C.H., and J.C. Wyngaard, (1989) "Evaluation of Turbulent Transport and
Dissipation Closures in SecondOrder Modeling." American Meteorological Society,
V.46, No.4, 2311 2330.
Nikjooy, M., H.C. Mongia and G.S. Samuelsen, (1989) "Calculations of Recirculating
Flows Using SecondMoment Closure." 27th Aerospace Sciences Meeting, January
912, 1989/Reno, Nevada, 11 p.
Oliver, H.R., (1971) "Wind Profiles in and above a Forest Canopy." Quarterly J. of
the Royal Meteorological Society., V.97, 548553.
Orszag, S.A., (1980) "Numerical Simulation of Shear Flows over Wavy Boundaries."
NASA Contractor Report 3271, 35 p.
Owen, P.R., (1964) "Saltation of Uniform Grains in Air." J. of Fluid Mech., V.20, Pt.2,
225242.
Perrier, E.R., J.M. Robertson, R.J. Millington and D.B. Peters, (1971) "Spatial and
Temporal Variation of Wind Above and Within a Soybean Canopy." Agricultural
Meteorology, V.10, 421 442.
Pinker, R.T., and J.Z. Holland, (1988) "Turbulence Structure of a Tropical Forest."
Bound aryLayer Meteorology, V.43, 4363.
Pinker, R.T., (1983) "The Canopy Coupling Index of a Tropical Forest." BoundaryLayer
Meteorology, V.26, 305311.
Raynor, G.S., (1971) "Wind and Temperature Structure in a Coniferous Forest and a
Contiguous Field." Forest Science, V.17, No.3, 351363.
Shaw, R.H. and I. Seginer, (1987) "Calculation of Velocity Skewness in Real and Artificial
Plant Canopies." BoundaryLayer Meteorology, V.39, No.4, 315332.
Shih, TH. and J.L. Lumley, (1986) "SecondOrder Modeling of NearWall Turbulence."
Physical Fluids, V.29, No.4, 971975.
Skupniewicz, C.E. R.F. Kamada and G.E. Schacher, (1989) "Turbulence Measurements
over Complex Terrain." BoundaryLayer Meteorology, V.48, 109128.
Spalart, P.R., (1988) "Theoretical and Numerical Study of a ThreeDimensional
Turbulent Boundary Layer." J. of Fluid Mechanics, V.205, 319340.
Speziale, C.G., (1989) "Discussion of Turbulence Modeling: Past and Future." Final
Report for Institute for Computer Applications in Science and Engineering, Report
No. 8958, 23 p.
4
I .
Speziale, C.G., (1982) "On Turbulent Secondary Flows in Pipes of Noncircular
CrossSection." International J. Engineering Sciences, V.20, No.7, 863872.
Thom, A.S., (1971) "Momentum Absorption by Vegetation." Quart. J. R. Met. Soc.,
V.97, 414428.
Thompson, J.F., Z.U.A. Warsi and C.W. Mastin, (1985) Numerical Grid Generation
Foundations and Applications Elsevier Science Publishing Co., Inc., New York, 483p.
Uchijima, Z. and J.L. Wright, (1964) "An Experimental Study of Air Flow in a Corn
PlantAir Layer." Bull. Nat. Inst. Agr. Sci., Ser. A, No.ll, 1965.
Wilson, J.D., (1988) "A SecondOrder Closure Model for Flow Through Vegetation."
BoundaryLayer Meteorology, V.42, 371 392.
Wilson, N.R. and R.H. Shaw, (1977) "A HigherOrder Closure Model for Canopy Flow."
J. Applied Meteorology, V.16, 11971205.
Wyngaard, J.C., (1987) "A Physical Mechanism for the Asymmetry in TopDown and
BottomUp Diffusion." American Meteorological Society, V.44, No.7, 10831087.
Wyngaard, J.C. and R.A. Brost, (1984) "TopDown and BottomUp Diffusion of a Scalar
in the Convective Boundary Layer." American Meteorological Society, V.41, No.l,
102112.
University of South Florida, Department of Geology and University of Florida,
Department of Coastal and Ocean Engineering, (1988) "Beach Monitoring Project
Sand Key Phase II Beach Nourishment Program (North Redington Beach and
Redington Shores PostNourishment" Profile Sediment Summary Data.
Zhang, J. and B. Lakshminarayana, (1989) "Computation of Three Dimensional
Turbulent Boundary Layers in Internal Flows, Including Turbomachinery Rotor
Blades." Proc. Ninth Intern. Symp. on Air Breathing Engines, V.1, 529538.
Appendix B
OneDimensional Canopy Flow Model
Appendix B
OneDimensional Canopy Flow Model
Introduction
The current version of the FDNR VEG code employs a onedimensional second order
turbulence model to determine the effects of the canopy on the flow and the shear stress at
the ground. Following is a review of the onedimensional model and its implementation in
the FDNR VEG program. Part A is a brief discussion of the turbulence model, numerical
solution procedure and calculation of the bottom shear stress. Part B is a discussion of the
model calibration, model parameters and model implementation.
Part A
Governing Equations
The onedimensional, second order closure model consists of five partial differential
equations for the five variables 5, u'w', u'u', v'v' and w'w'. These equations have been
developed with the assumption that the flow is steady and the canopy is uniform in the
horizontal direction. The equations are:
Of 9 u'w'
1) + + + CdAU = 0,
at az
8u'w' 9 w q
2) + w'w' 2 + u'w Cq2 = 0,
9t + z Qz Qz 32 Qz
9u'u'   9 9 u' q ( q 2q3
3) + 2u'w' qAz ) + u'u'2 + 3 = 0, (1)
) v'v' a v'v'\ q ( 2q3
4) a qA V'V' + = 0,
at 9z 8z 3A2 3 3A3
5) w'w' ( z ww' q (2 2q3
5) qA1 )+ w'w' + = 0,
at az 0z 3A2 3 3A3
where z is the height, t is time, A1, A2 and A3 are length scales, C is a constant, Cd is the
canopy drag coefficient, A is the leaf area density, and q is the square root of the turbulent
kinetic energy per unit mass,
q = (u'u' + 7vv + w'w')1/2
The three length scales AX, A2 and A3 are defined in terms of a single length scale, :
A1 = 81i, A2 = 862 and A3 = 63,
where 86, 62 and 83 are constants. The length scale I is defined as kz but is then modified
due to the presence of the canopy using the approach outlined by Wilson and Shaw,
(1977). The drag coefficient Cd and leaf area density A are properties of the canopy, the
constants C, 81, 82 and 83 are model parameters.
The canopy parameters Cd, A and the model parameter 81 are discussed in Part B. The
remaining model parameters, C, 82 and 83, can be obtained by requiring the model to
exactly reproduce a constant stress layer in the absence of a canopy.
Empirical Constants
The method to determine the constants C, 82 and 63 is as follows: For a constant stress
layer, the following conditions must be met,
u'w' = u2 and
9u u.
= (2)
,z Kz
where u. is the friction velocity and t is von Karmon's constant. Also, the following
relationships have been determined from measurements in a constant stress layer,
u'u' = 3.5 uw ,
v'v' = 1.5 u'w' and (3)
w'w' = 1.5 w'w' .
Upon substituting Equations 2 and 3 into Equation 1, three linear algebraic equations for
the three constants can be obtained. Their solution yields the following values for these
constants:
62 = 0.85,
63 = 16.57 and
C = 0.077.
Boundary Conditions
The governing equations are valid in a region from the ground to some height above the
canopy, where a constant shear layer is assumed to exist. In order to solve these equations
in this region, boundary conditions must be specified. At some distance above the canopy,
H, the following conditions are enforced:
u = UH, = 0,
= 0, =0 and
9z 9z
w'wt'
0.
9z
At the ground, the velocity is assumed to be zero at a height zo (i.e. at the height of the
bottom roughness elements). Here, the conditions
U = 0,
9u'w' 9u'u'
a = 0 a O,0
az &z
vv'1 Ow'w'
= 0, and 0,
9z 9z
are enforced. Implicit in these conditions is that a constant shear stress layer exists near
the bottom boundary.
The value of zo used to define the bottom roughness depends upon the flow conditions,
specifically the value of the shear at the bottom. Owens (1964), has developed a
relationship between zo, the critical stress for motion and the bottom stress (friction
velocity).
SZn U* < U*c
zo u2 (4)
2U U*>U*C
where # is an empirically determined coefficient, g is the acceleration due to gravity and z,
is the bottom roughness based on the grain size. This formulation, which is used in the
onedimensional model, alters the bottom roughness in the presence of a saltation layer and
is a generalization of the results of Bagnold (1941).
Numerical Solution Procedure
An implicit second order accurate numerical solution scheme is employed to solve the
governing equations. For the purpose of developing the finite difference scheme, the
equations are rewritten in vector form as:
oQ, aF
+ + = 0 (5)
9t 8z
where
UWW
Suu'w'
u W q A 9 u w
S uu' F \ u'u'
v'v"' 8* O~v
W'W'l
3qA w'w'
and
CdAU2
ww'u u'w' + Cq2a
2) 2 2
,= 27w'  (u'u' q)
2) 2 2
2(w'w q 2 
Equation 5 represents a system of five equations for five dependent variables contained in
the vector Q. To obtain a numerical solution, the vertical domain is represented by
discrete points, zi, which span from zo to H. Between each grid point the spacing Az is
defined as
Azi+1/2 = Zi+1 zi.
In order to provide adequate resolution of the near ground shear layer, the grid points were
clustered near the bottom using a hyperbolic clustering function (see Thompson et al.,
1985). In general, the spacing between the first and second grid point is of the order (zo).
The elements of Q are then approximated by discrete values, Q, at each grid point i.
Using this approximate representation, the spatial derivative of the term F is
approximated with finite differences:
SF F F
S F (1+1/2) (1/2)
z A/z
where
F + F
F (1+1)
(1+1/2) 2
Derivatives of the variables in F and S are approximated in a similar way.
Using these finite difference approximations, Equation 5 is transformed into a nonlinear
coupled system of first order differential equations for the 5 N variables, Q :
9Q
+ 0 + =o0. (6)
at
The steady state solution to Equation 6 is the asymtotic solution to the unsteady
equations. To obtain this asymtotic solution, an initial guess is made for the 9 for each i.
Equation 6 is then marched in time until it converges. An implicit time marching scheme is
used, i.e.
AQ + F(+) S(n+i) =0, (7)
At
where
i i i
and n indicates the time step.
Equations 7 are nonlinear due to the terms (n+ and n+) (they require values for
i i
(n+ which are yet unknown). Therefore, the terms F and S must be linearized in
time:
a F
(n+i) S" S
aF aS
where  and are the Jacobian matrices, herein denoted as G and H, respectively.
i i
I
Substituting the linearizations of and into Equation 7, we get
i i
A + (GA ) +HAQ =6 F+
At i + HA Q i
which can be rewritten as
[I + At (,Gi + Hi)]A Q =At (, F + S (8)
This equation is a matrix system of equations for the dependent variables A Q at the grid
point z;. Here I is a 5 by 5 identity matrix. Note that the right hand side of Equation 8
represents the steady state solution. When the right side approaches zero, A Q + 0 and
the steady state solution is obtained.
Equation 8 can be written for each grid point i, producing a system of 5(N 2) equations.
These equations are coupled, however, since each set of 5 equations at the ith location (i)
are related to neighboring equations at the (i + 1) and (i 1) positions through the spatial
derivatives. As a result of this coupling, the system of equations has the form of a block
tridiagonal system with 5 x 5 blocks. These equations can be solved efficiently using a
block tridiagonal solver analogous to the Thomas Algorithm developed for scalar
tridiagonal systems (see Anderson et al., 1984).
After solving the system of equations the dependent variables for all grid points are
updated:
(n+l) n
This procedure, which constitutes one complete time step, is repeated until the solution for
9 converges.
The steady state solution provides vertical profiles of the mean flow i and the four
2 2 .2
Reynold's stresses u'w', u' v' and w2 The shear stress at the ground, which is required
to determine the sand transport, can be obtained from the steady state solution as follows:
rg = pu.
where u. = (u'w')z.
Here 9 is the ground shear stress, p is the air density, u. is evaluated at the ground (i.e. at
z0). The ground shear stress rg is then used to determine the sand flux.
j
Part B
Model Parameters
The onedimensional canopy flow model requires the specification of three parameters.
Two of the parameters, the plant drag coefficient (Cd) and the vertical leaf area density
profile (LAD), characterize the plant canopy and its effect on the flow. The third
parameter, 81, is a coefficient which scales the turbulent transport of the Reynold's stresses.
Each of these parameters is expected to vary with the type of plant canopy. In addition,
the magnitude of the LAD will vary as the plant spacing and height are varied.
The LAD and the Cd are usually obtained from measurements or experiments on a
particular canopy, whereas the coefficient 81 is obtained by calibrating the model with
experimental mean wind velocity and LAD profile data. The available experimental data is
adequate to estimate the plant LAD profile but is insufficient to determine the drag
coefficient. Therefore, Cd was also obtained from model calibration. The methods used to
determine the three parameters LAD, Cd and S1 for the plant types contained in the
database are discussed below.
Leaf Area Density
The leaf area index (LAI) is a measure of the onesided leaf area of a plant divided by the
total ground area surrounding the plant. The LAD is the LAI per unit height, which can
be thought of as the vertical distribution of leaf area with height of the plant. At present,
there are no direct measurements of the LAD for Florida's coastal plants. However, these
values have been either measured directly or estimated for various common cultivated
plants. Therefore, the LAD for coastal plants has been estimated based on the data
available for noncoastal plants.
There are a number of published studies of turbulent flow in canopies for crops, orchards
and forests due to its importance in the understanding of photosynthesis, evaporation, heat
transfer and aerosal dispersal within these canopies. Ten such studies have been found in
the literature that include measurements of the vertical leaf area density profile. These ten
plants are used in this study as "reference" plants and the coastal plants were "matched"
to these plants. The plants' physical characteristics, which were obtained from botanical
literature, botanists, and photographs, were used as a basis for comparison. Specifically,
such features as the general shape of the plant, the leaf texture (size, shape and number of
leaves), and the overall flexibility of the plant were compared. The reference plants and
their physical characteristics are listed in Table B1 and the corresponding leaf area density
profiles are shown in Figures B1B10. Table B2 lists the coastal dune plants under
consideration and their physical characteristics.
General application of the FDNR model requires values for the LAD for any plant spacing
or height. Unfortunately, the published data usually contains data for one set of conditions
and information, i.e. one spacing and one height. In order to obtain LAD profiles for plant
spacings and heights, other than for the reference condition, the following expression was
used:
LAD(z, L, H) = LAD(z, Lo, Ho) (H/Ho) (Lo/L)2, (9)
where Lo and Ho are the reference plant spacing and height, and L and H the actual plant
spacing and height, respectively. The H/Ho term proportions the actual plant height to the
reference plant height. So, if the spacing remains the same, but H is smaller than Ho, then
the LAD decreases. The Lo/L term scales the magnitude of the LAD profile depending on
the actual spacing relative to the reference spacing. If the height remains the same but the
spacing increases, the LAD(z, L, H) decreases because of the increased area, L2, associated
with the plant.
Drag Coefficient, Cd and 61
As stated previously, the drag coefficient is usually obtained from measurements on a
particular plant, and 61 is calibrated from experimental data. Because of insufficient
information in the literature on Cd, both Cd and 61 had to be determined so as to produce
the best fit between computed and measured data for u.
The Cd and 61 for the reference plants were calibrated with the onedimensional model and
the measured wind velocity profiles obtained in the literature for the reference plants. The
Cd and 61 were varied in the model to produce velocity profiles that best fit the existing
measured wind velocity profiles. The parameter 61 was given just two values, 1.0 and 0.1,
while the Cd was changed to produce the best fitting profile for each S1. The model profiles
compared with the experimental profiles for each reference plant are shown in Figures
B11B20.
Table Bl:
REFERENCE PLANTS
Plant Leaf
Vegetation Type Height Length Spacing LAI Cd Flexibility
(cm) (cm) (cm)
Almond tree 762 8.9 800.0 1.3 .15 low
Broad Bean shrub 118 8.9 76.2 6.25 .07 mod
Corn grass 140 61.0 40.0 4.2 .17 mod
Oak tree 2300 15.2 243.8 4.9 NA low
Sweet Orange tree 425 10.2 400.0 6.5 NA low
Pine tree 1050 8.9 260.0 NA NA low
Rice grass 76 22.9 21.0 3.8 0.210 high
Soybean shrub 108 11.4 76.0 4.5 0.035 low
Sitka Spruce tree 1050 1.9 262.0 10.0 0.140 low
Common Wheat grass 125 29.2 18.0 3.0 0.200 high
NA = not available
Table B2:
COASTAL PLANTS
Leaf
Vegetation Plant Height Length Spacing LAI Flexibility
Type (cm) (cm) (cm)
Dune Prairie Grass grass 150 137.0 46 4.60 high
Saltmeadow Cordgrass grass 90 61.0 12 5.00 high
Sea Oats grass 183 60.0 60 4.00 high
Beach Bean herb 15 7.6 122 3.84 mod
Beach Croton herb 70 5.2 60 1.68 mod
Beach Morning Glory herb 15 10.0 27 4.32 mod
Dune Sunflower herb 90 7.6 76 1.60 mod
Railroad Vine herb 15 10.0 21 4.90 mod
Sea Blite herb 76 5.2 61 2.40 low
Sea Purslane herb 15 3.8 30 1.20 mod
Camphorweed shrub 70 5.2 90 1.25 mod
Cocoplum shrub 300 5.1 460 10.00 low
Conradina shrub 90 1.2 150 1.20 mod
Gray Nicker Bean shrub 150 5.2 61 7.10 mod
Prickly Pear Cactus shrub 90 21.0 30 2.40 low
Rosemary shrub 90 1.2 180 1.20 mod
Saw Palmetto shrub 150 76.0 15 4.20 mod
Sea Grape shrub 240 20.0 240 10.00 low
Sea Myrtle shrub 200 3.8 150 3.60 mod
Spanish Bayonet shrub 215 52.0 300 3.20 low
Australian Pine tree 3000 21.0 600 23.60 low
Bay Cedar tree 2400 3.8 150 10.60 low
Cabbage Palm tree 2300 75.0 600 7.10 low
Chapman's Oak tree 600 7.6 180 15.10 low
Coconut Palm tree 900 460.0 750 12.60 mod
Myrtle Oak tree 300 7.6 120 9.80 low
Sand Live Oak tree 900 5.2 460 18.90 low
Sand Pine tree 900 7.6 460 12.60 low
Slash Pine tree 1830 25.0 900 18.80 low
Southern Magnolia tree 180 15.0 900 14.10 low
Southern Red Cedar tree 760 6.0 600 15.70 low
Wax Myrtle tree 520 10.0 280 3.50 low
Leaf Area Density Profiles of Reference Plants
Almond
H. = 26.3 ft.
L. = 26 ft.
Figure B1
Corn
H. = 4.6 ft
L. = 1.3 ft
1.20
1.00
0.80
20.60
0.40
0.20
Oak
H. = 75.5 ft
L, = 8.0 ft
Leaf Area Index
Figure B4
1.20
1.00
0.80
S0.60
0.40
0.20
1.20
1.00
0.80
20.60
0.40
0.20
Figure B2
1.20
1.00
0.80
o
=0.60
0.40
0.20
Figure B3
Leaf Area Density Profiles of Reference Plants
Orange
H. = 13.9 ft
. = 13.1 ft
1.20
1.00
0.80
30.60
0.40
0.20
Pine
H. = 34.5 ft
L. = 8.6 ft
Leaf Area Index
Figure B5
Rice
H. = 2.5 ft
I.. = 0.7 ft
Figure B7
Leaf Area Index
Figure B6
Soybean
1.20 H. = 3.5 ft
= 2.5 ft
0.80
I 0.60
0.40
0.20
0.00 ..... ....... ......................... o
0.00 2.00 4.00 6.00 8.00 10.00
Leaf Area Index
Figure B8
1.20
1.00
0.80
0.60
0.40
0.20
1.20
1.00
0.80
=I0.60
r
0.40
0.20
Leaf Area Density Profiles of Reference Plants
Spruce
1.2o H. = 34.5 ft
L, = 8.6 ft
1.00
0.80
IO.60
0.40
0.20
0.00 .. .. .
0.a0 0.10"o.o 0.oo0 0.4
Leaf Area Index
Figure B9
Wheat
1.2o H. = 4.1 ft
L, = 0.6 ft
1.00
0.80
= 0.60
I c.o
I
0.40
0.20
0.00 ..... .
0.00 4.00 8.00 12.00 16.00
Leaf Area Index
Figure B10
I
Experimental and Model Wind Velocity Profiles
of Reference Plants
Almond
+ / /
dd = 1.0
2.00 4.60 6.00
Mean Velocity
Sr
2.00
1.50
0.50
2.00
1.50
2 1.oo
0.50
0.00
0.00
Bean
I
+
+++++ experimental
d, = 0.1, C,
d, = 1.0, Cd
4.
nT1n
0.20 0.40 0.60 0.80
Mean Velocity
1.00 1.,
Figure B11
Figure B12
Corr
/
+f
+++++ experimental
 d = 0.1, Cd
S  d, = 1.0, Cd
I 0.20
.m
data
= 0.15
= 0.17
2.00
1.50 0
1.50
r
o1.oo
"I
0.00
0.01
Oak
I
I
S+
+++++ experimel
+ d = 0.1
d = 1.0
, ,, ,, vn ............. I In....I
2.00 4.00 6.00 .0 1.0
Mean Velocity
Figure B13
ntal data
, Cd = 0.18
, C = 0.29
8.00
a)
0
'5
U.U0 +r
0.00
data
= 0.06
= 0.065
io
2.00
1.50
?1.00o
0.50
0.00 
0.00
i... 1[1... 11. 111. 11 =l ll, T, ....n ,. ,
0.40 0.60 0.80 1.00 1.20
Mean Velocity
ntal data
SCd = 0.05
. C = 0.06
Figure B14
l i...... i lll.. i. li.l.. .. l.l. ll.. l
,
I
1
0
Experimental
and Model Wind Velocity Profiles
of Reference Plants
Oronae
/
/
/
/
//
+ +++++ experim
d = 0.
d11.
I
2.00
1.50
r
.?1.0o
0.50
2.00
1.50
.1.00
0.50
0.00
0.00 o ,,,, ,, l,, ,,, ,,,1i. 0,,,,,,
0.o 00 6:" .o i: 2o 2
Mean Velocity
Figure B15
Rice
+
+ /
S +++++ experimental data
S d,= 0.1. C, = 0.14
S d, = 1.0, Cd = 0.18
. ...... '... .. '.2...... 0.o
0.50 1.0Velci
Mean Velocity
2.00
1.50
1.00
0.50
0.00oT
0.00
Pine
+++++ experimental data
d d = 0.1, Cd = 0.06
d, = 1.0, Cd = 0.075
Mean Velocity
Figure B16
Soybean
+++
Sepeiena dt
experimental data
d, = 0.1, Cd = 0.10
d, = 1.0, Cd = 0.19
200 4.00 6. 8.00
Mean Velocity
Figure B17
mental data
1, C, = 0.6
0, Cd = 0.5
.0
5
2.00 
1.50
1..oo
0.50
.4
a
.I
I
0.00o
A 
 '
*t
Figure B18
Experimental and Model Wind Velocity Profiles
of Reference Plants
Spruce
+ 
+
+++++ experime
S d = 0.1
S d, = 1.0
I / *
'1 16 6 .. Y.. 6. . ..
1.00 2.00 3
Mean Velocity
.00b
ntal data
. Cd = 0.5
i, Cd = 0.7
... 0
4.00
2.00
1.50
r
.21.00
0.
0.50
Wheat
/
/
/
+++++ experimental data
d = 0.1, Cd = 0.09
d, = 1.0. C = 0.12
Mean Velocity
Figure B20
2.00
1.50
01.oo0
0.50
0.
0
I
0.00
0.00
Figure B19
Appendix C
TwoDimensional Canopy Flow Model
Appendix C
TwoDimensional Canopy Flow Model
Introduction
A twodimensional second order turbulence model for flow through a finite length canopy
on a flat terrain has been developed in order to determine twodimensional effects such as
the dependence of the transition length on the canopy properties. Following is a
description of the twodimensional model including the governing equations, boundary
conditions, and numerical method used to solve the governing equations. Intermediate as
well as final (parameterized) results are presented and discussed.
Governing Equations
The twodimensional second order turbulence model consists of six partial differential
equations for the six dependent variables u, v, u'w', u'u', v'v', and w'w'. They are:
S+ 0 = 0, (1)
ax Oz
u O9u2 9 w Ou 6 u'w' aP
+ + + + = CdA UIl, (2)
at xz 6z ox dz ax
9W aW OW2 u'w w'w' W P
+ + + + +  + = CdA TI, (3)
Ot + 9x Q+ W z
+ a )+u 2) W o
+ iW + W Cq2) 8u + ('' Cq2) + U' W'
_9z az ax
a au'w' a aa0 9w'w'
OqAz u' 2 z q qAOx qAz
9x ax az az ) ax z az ax
Su'w', (4)
3A2
u'u' au'u' au'u' 2 u 
5t +  +  + 2 (uu Cq2) + 2uw
at x az x Oz
3 u g qA 2 qA
ax ax Qz z ) az 'axj
=3 3 (5)
312 \ 3d3'
dvW __OvW' _9v'v' 9 / 9vv"\ 9 I O v'"
_a7 P a V d a (qX ____
at + ii + U7  ql
xt a z ax ax az az9
3 v'3 and (6)
3A2 3 3A3
aw'w' aw'w' _aw'w' a 
+ + w +2u'w' +2 w'w' Cq2) 
q i 3 q __ q3
q / q2v 2q 3
S3A 3 3A3'
where z is the height above the ground and x the horizontal downstream axis. The
turbulent kinetic energy/unit mass, q2, and the parameters C, Cd, 61, 82 and 63 are defined
and discussed in Appendix A. These equations are similar to those of the onedimensional
turbulence model except that: 1) the equations for the Reynolds stress have advective and
diffusion terms in twodimensions, 2) there is an additional momentum equation for T (the
vertical velocity component), 3) the continuity equation has been included, and 4) both
momentum equations include the dynamic pressure P.
Boundary Conditions
Boundary conditions for the governing equations must be specified at the edges of the
domain, that is at the ground, at some height above the canopy, at some distance upstream
of the canopy and at some downstream point. The top and bottom boundary conditions
are identical to those in the onedimensional model with the addition of conditions for w
(zero at both the ground and the upper boundary). At the upstream boundary, all the
variables are specified for a bare ground constant stress layer. At the downstream
boundary, all variables are extrapolated from the adjacent interior domain. Note that the
downstream boundary can be placed in the canopy (for very long canopies) or on bare
ground beyond the canopy (for shorter finite length canopies). Implementation of the
boundary conditions, however, is not affected by this choice.
Numerical Solution Procedure
The method employed to solve the twodimensional second order turbulence model is a
semiimplicit finite difference scheme. It is however, more complicated than the scheme
employed in the onedimensional model (Appendix A) due to two factors. First, since the
equations are twodimensional, a spatially factored scheme is necessary to efficiently solve
the Reynolds stress and momentum equations. Secondly, since the momentum equations
I
contain the dynamic pressure term and the continuity equation is not automatically
satisfied (as it is in the onedimensional model), a special time splitting technique (see Kim
and Moin, 1985) and staggered grid is applied to solve the momentum equations.
Described below are: 1) the grid system, 2) the spatially factored scheme for the Reynolds
stress equations, and 3) the method used to solve the momentum and continuity equations.
Grid System
The twodimensional grid system employed to represent the domain is a rectalinear system
with variable spacing in one direction and equal spacing in the other direction. The
constraint of equal spacing in one direction is due to the solution method of the momentum
equations and is discussed in a later section. Variable spacing was used in the vertical
direction in order to cluster grid points near the bottom and resolve the near ground shear
layer. Figure C1 shows the grid system and variable location.
Note that the Reynolds stresses are defined at the nodes, the mean velocity
ij
components are "staggered" with the uij and vij components defined on the vertical and
horizontal cell faces respectively and the pressure Pij is defined at the cell centers. Also
defined on the nodes are values of the leaf area density Aij. The i and j subscripts refer to
the x and z grid points respectively.
Numerical Solution of the Reynolds Stress Equations
Development of the spatially factored scheme proceeds first by rewriting the four equations
(Equations 47) for the stress components in vector form:
aR aF aG
o 4 (8)
at 9z Oz "'
where
"'U' 3q '()
F (10)
uv' v qA q'
ww q 2q1 u'w'
\ L/J _FZ_ I
i
wu'w 2q\A X 1 q 'W "
W7i~s qAj 6FU 2gA^' UI
~= v and (11)
ww'w' 3qAiX
u 1 w rw, _uw
Y(wO aw Cq()T w W'
22) U . (u'u' 2
2u'w'w ( 2 (w'w Cq2) (ww 
An implicit time marching algorithm (in delta form) is first introduced and can be
represented as:
3 + 5 (AARij) + Sy BA R) R A R = (RHS)1, (13)
At
where ARij = R+1 Rin. (n indicates time level), A and B are the Jacobian matricies
8F oB
A = and B = 
a R QR
andRHSi, is n + .
13 S3 13
The difference operators 86, 6y and other derivatives are the same as those used in
Appendix A. Note that the fluxes at the n + 1 time level have been linearized in time, i.e.
SF
e = r +F a A R. (14)
Equation 13, when written for each grid point i, j represents a coupled system of
equations, forming a sparce banded matix which can be solved at each time step to obtain
the change in the Reynolds stresses at each point A R The new values of the stresses
ij
can then be obtained from
R + Rn + R (15)
ij ij ij
The large (IMAX X IMAX x 4) sparse banded matrix must be inverted each time step and
solved approximately using a spatially factored scheme. This approach is much more
efficient than a direct solution of the system. The spatially factored scheme is expressed as
I + AtAi At (I+ AtSBi) A = At R (RHS),. (16)
Note that the difference between this equation (when expanded) and Equation 13 is the
a S
term At2(5xA )SyB. This term is of order At2, thus only a small error is introduced
SR
by omitting it. Also, once the converged solution is obtained, ARij = 0, RHSy = 0 and
subsequently the error will also go to zero. Thus this error will not affect the steadystate
solution.
Equation 16 can be written as a two step process:
(I + AtSA) AR!iAtRHSij and (17)
(I + AtsyBi,) ARi = ARjj. (18)
This splitting of the operator requires simple tridiagonal matricies to be solved for all i
lines and then for all j lines.
The solution procedure for the Reynolds stresses consists of the following steps.
1. Assume an initial Reynolds stress profile (usually a constant stress layer).
2. Form Equations (4) and solve for ARij (using current values for U and u).
Dn+1
3. Obtain n+ for A R
ij i3
4. Update the boundary conditions.
5. Repeat steps 23 until ARij + 0.
Note that the term S contains values of the mean velocity components. These values
t i
come from the solution of the momentum equations described below.
Numerical Solution of the Momentum Equations
The solution procedure for the mean momentum equations and continuity (Equations 13),
requires special treatment due to the presence of the dynamic pressure term. The scheme
used here is similar to that of Kim and Moin (1985) and is briefly described below:
The governing momentum equations can be written as:
9u _9u2 8uw 9u'u' 9u'w',
T + P =  CdAl` uiu and
9t 9x 9z 9x 8z
(19)
U_ aUV V2 auw' w'w' U.
t x 9z Ox oz
First by defining the following intermediate velocities W and U* as
9u* 9
= + P, and
(20)
Ot Ot
and then applying standard finite difference approximations to Equations 19, 8t and 
are determined at each grid cell face.
Next, the pressure field is determined such that the velocity field is a divergence field. The
equation for the pressure is obtained by substituting Equation 20 into the continuity
equation:
(a O i ifO= a 8W \ _(2PO2P (21)
Ot' 9x + z) =9t x + z 9x2 + z2J (21)
By requiring that the velocity field remain divergence free, (left hand side of Equation 21 =
0) an elliptical partial differential equation is obtained which can be solved for the pressure
field. When solved, the pressure field can be substituted into Equation 20 and and
obtained, which are then used to advance the velocities in time.
The numerical procedures used here are described by Kim and Moin, (1985). Note that the
elliptical equation for pressure is solved directly using transform methods which requires an
equally spaced uniform grid in one direction. We have modified the approach to employ
Fast Sine and Cosine Transforms so that the solution procedure is efficient. This
modification further restricts the number of grid points used to be 2N + 1 where N is an
integer.
The momentum equations are solved each time step along with the Reynolds stresses and
the time marching scheme is continued until a steady state solution for the Reynolds
stresses and mean velocity components are obtained.
Sample Results
The twodimensional second order turbulence closure scheme has been applied to a finite
length canopy in order to demonstrate its use. The sample case consists of a 2 m high,
64 m long canopy. Its dimensions and coordinate system are shown in Figure C2. The
canopy has a triangular vertical LAD profile (Figure C3) with an LAI of 0.1 and is uniform
in the horizontal direction.
The resulting velocity profiles calculated using the twodimensional model are shown in
Figure C4. Upstream of the canopy a loglayer profile exists over the bare ground (Figure
C4a). At the beginning of the canopy the flow is retarded at peak value of the LAD
vertical profile but jetting of the flow occurs due to the lower LAD values near the ground
(Figure C4b). Further into the canopy the jetting subsides. The near ground velocity and
stress continue to decrease until the transition length is exceeded. At this point (Figure
C4c) the velocity profile is the same as that which would be produced by the
onedimensional second order turbulence closure model and the bottom stress reaches its
minimum value. This profile exists until the end of the canopy is reached (assuming, of
course, that the canopy is longer than the transition length). Downstream of the canopy a
wake region exists (Figure C4d) over the length of which the velocity profile returns to the
upstream loglayer profile. In this region the bottom stress continually increases until the
bare ground value is obtained.
TwoDimensional Effects
As seen in the sample case above, a number of twodimensional effects can occur which will
affect the total transport of sand into and out of a canopy. The twodimensional model has
been used to predict the twodimensinal effects that occur in finite length canopies and the
results have been incorporated in a second version of the FDNR VEG code (Version 2.0).
There are three lengths which have been defined to characterize the twodimensional
effects. The lengths are depicted in Figure C5 and are defined below:
Li: The total extent of the disturbed air flow measured from the beginning of the canopy
to the point in the wake where the upstream velocity profile is restored.
L2: The distance into the canopy over which sand is deposited.
L3: The length in the wake region over which sand is eroded.
The value of these lengths depends on the LAI of the canopy and the length of the canopy.
The value of L2 and L3 also depend on the value of the critical shear stress relative to both
that of the bare ground and the minimum value on the bottom in the canopy.
Note that the minimum value of the bottom stress in the canopy is itself dependent on the
length of the canopy. In cases where the canopy is longer than the transition length, the
bottom stress decreases to a minimum value which is equivalent to the value obtained using
the onedimensional model. However, if the canopy is shorter than the transition length
then the minimum bottom stress will not decrease to the onedimensional model value.
The length L1 includes both the canopy length and the length of the wake region. Results
from the twodimensional model indicate that the wake region will range from 22 to 45
canopy heights, thus L1 will range from 22 to 45 canopy heights plus the length of the
canopy. In general, the wake region increases with decreasing LAI and increases with
increasing canopy length until the canopy length exceeds the transition length. Beyond this
point, the wake region length remains constant.
The values of L2 and L3 have a similar dependence on LAI and canopy length, but also
depend on the value of the bottom stress. Two distinct trends for L2 and L3 are apparent
depending on whether or not the canopy reduces the bottom stress to a value below the
critical stress. If the canopy does not reduce the bottom stress to a value below the critical
stress, then sand is transported throughout the length of the canopy. In this case, L2 is
simply the length of the canopy and L3 is the length from the end of the canopy to the end
of the wake region, i.e. L3Lc, where Lc is the canopy length. The net rate of sand
accumulation in the canopy is then the transport rate upstream of the canopy (i.e. bare
ground rate) minus the transport rate at the end of the canopy. If the canopy reduces the
bottom stress to a value below the critical stress, then sand transport stops somewhere
within the transition length of the canopy and does not intiate again until the bottom
stress rises above the critical stress in the wake region. In this case L2 is less than Lc and
L3 is less than the wake region length. In general L2 and L3 decrease with increasing LAI
and increase with increasing Lc until Lc surpasses the transition length. At this point L2
and L3 are constant.
The magnitude of the wind may also affect the value of L2 and L3. At lower wind speeds,
the canopy can reduce the bottom stress to below the critical value sooner than for higher
wind speeds. For certain canopy LAIs and lengths, changes in wind speed may cause
changes in the values of L2 and L3. However, the length L3 is essentially independent of
wind speed.
The data obtained using the twodimensional model has been incorporated into a database
for use in the FDNR VEG (Version 2.0) code. Results from the onedimsensional model,
the canopy length, height and the width are used in conjunction with the database to
determine the minimum value of the bottom stress that occurs. Subsequently the
characteristic lengths L1, L2, L3 and the net sand accumulation rate over the area of the
canopy are calculated. Note that use of the onedimensional model results include such
effects as bottom roughness and variations in LAD profiles, thus they do not have to be
included as part of the twodimensional database.
dx
i +1 **
xPj dz
j+1
Figure C1 Staggard grid system and variable locations.
4.00
N
0.00
40 20 0 20 40 60
x/H
Figure C2 Canopy dimensions and coordinate system.
I
2.0
1.5
4
I
1.0
0.5
0.0
0.00 0.05 0.10 0.15
LAD
Figure C3 Canopy leaf area density profile (LAD).
Top of Canopy

5
Velocity
10 15
(at x/H = 20)
3
Top of Canopy
2 
0 5 10 15 2(
Velocity (at x/H = 0.5)
S2 Top of Canopy
N2 
Velocity (at x/H =
Top of Canopy
~ ~ I    
5
Velocity
10 15 20
(at x/H = 45)
Figure C4 Velocity profiles for 2D flow through canopy.
deposition
deposition
 iL2
deposition
IL2
erosion
La 
Figure C5 Definition of lengths L1,L2,L3 and Lc.
_____ ^
Appendix D
Flow Over an Irregular Boundary
Appendix D
Flow Over an Irregular Boundary
The first step in the development of a twodimensional model for the wind flow and sand
transport through a vegetation canopy on an irregular (wavy) terrain has been completed.
This flow situation is a generalization of the twodimensional flow through canopies on flat
terrain. The flat terrain problem was addressed using the twodimensional model developed
in Appendix C. The approach taken here to develop a twodimensional model for an
irregular terrain was to start with the flat terrain model and make the nessessary
modifications to extend its applications to more general surface geometries.
This extension to general geometries is a complex task due to the curvature of the
boundaries. The primary difficulty lies in treating the boundary conditions, since for a
normal orthogonal rectilinear grid network, the grid points will not coincide with the
curved boundary surface. The standard approach to circumvent this problem is to develop
a boundary fitted coordinate system and transform the governing equations (Equations
C1C7) onto this new coordinate system. The generation of the grid system itself requires
the use of computation techniques but is a straight procedure for smooth geometries (see
Thompson, 1985).
The transformation of the governing equations onto the general curvilinear coordinate
system, however, is a tedious task and greatly increases the number of terms in the
governing equations (see Anderson, 1983). Additional terms appear due to the curvature of
the system and the nonorthonality of the coordinates. As a first step in developing the
twodimensional canopy flow model for general surface geometery, the computational
solution procedure has been developed for the momentum and continuity equations. This
initial model represents viscous laminar flow over an irregular boundary.
The governing equations for this flow are the viscous, incompressible NavierStokes
equations:
a = 0, (1)
9u a9u2 aw QP 8( 5 992
7+ + + = v ( + 2) and (2)
axm awa2 P2 a2
t + + O + z = Vo2 + (3)
at ox z z 9x2 9z2
These equations are similar to equations C1C3 but with the Reynolds stress terms
replaced with viscous stress terms. In developing a solution procedure, these equations
L
have been transformed onto a general curvaliner coordinate system and solved using finite
volume methods. The numerical algorithm is essentially similar to the timesplit scheme
used in the flat terrain model. However, since a curvilinear coordinate system is used here,
the direct solution procedure for the pressure Poisson equation is not applicable. Therefore,
a standard ADI line iterative scheme is currently being used. It should be noted that the
solution of the pressure Poisson equation, which is required each time step, requires a
significant amount of CPU time. It is therefore suggested that a more efficient multigrid
approach be implemented (see Briggs, 1987) when a turbulence model is introduced.
The viscous flow model has been coded and applied to flow over a wavey bottom to test
the algorithm. Figure D1 shows the resulting velocity field calculated by the model. In
Figure D4a, it is apparent that flow separation is imminent. In Figure D4b, corresponding
to a higher Reynolds number, separation has occurred in the trough region. These results
indicate that the model is functioning properly and will provide a good basis for extending
the model to include turbulent flow and canopies. It is however, important to improve the
efficiency of the pressure Poisson equation solver and to further test the model before
proceeding.
(a) No separation
(b) Flow separation
Figure D1 Velocity profiles for flow over a wavy bottom.
Appendix E
Case Study
Appendix E
Case Study
The following example illustrates how the VEG 1.5 and VEG 2.0 programs can be used to
determine the relative effects of different plant spacings, types and heights on sand erosion
and deposition within and behind existing and proposed coastal vegetation.
A developer in Redington Beach, Pinellas County has just purchased a parcel of west
facing beachfront property, 300 feet shorenormal by 1500 feet shoreparallel. The property
is bounded to the east by a 12 feet wide asphalt drive, and the area within and on either
side is dominated by 80 feet tall Australian pines spaced 15 feet apart (see Figure El). The
developer, a conservationist, would like to remove the exotic plant species and replace them
with native coastal vegetation, namely, sand live oak and dune prairie grass. The developer
wants to plant 10 feet tall sand live oaks with a 15 feet spacing and 0.25 feet tall dune
prairie grass 4 feet apart.
The main concern is that replacing the established plants with native vegetation may cause
sand to erode at a greater rate. The FDNR VEG program can be used here to compare the
sand transport through the nonnative planting with the sand transport through the
proposed native planting to determine the effects of the new planting on the sand accretion
rate. The program can also be used to determine, if necessary, what plant spacing and/or
height of the desired native plants will be necessary to produce the same sand trapping and
retention characteristics as the existing vegetation.
The onedimensional code, VEG 1.5 assumes that the flow is fully developed over most of
the canopy. It computes the rate of sand transport (in yds3/yd/hr) over bare ground and
within the canopy. The rate of sand accumulation within the canopy is, thus, the difference
in the bare ground and the canopy transport rates multiplied by the width of the canopy.
The twodimensional code, VEG 2.0, computes the rate of accumulation directly and
accounts for the transition regions.
Note that the canopy depth (300 feet) to height (80 feet) ratio in this example is such that
a fully developed flow will probably not occur within the existing canopy. Thus, VEG 1.5
cannot be used realistically for this case. Even though the proposed plant height is much
shorter (10 feet), it too may be borderline for its use in VEG 1.5. To illustrate these facts,
both versions (1.5 and 2.0) have been used to compute the rates of sand accretion within
both the existing and proposed canopies. The results are shown in Table El.
In this example, the grain size diameter is 0.181 mm and the sand gradation is uniformly
graded (USF, 1988). The wind speed for bare ground sand transport is 12.7 mph. Local
climatological data gathered from NOAA (National Oceanic and Atmospheric
Administration, Tampa, FL, 1988) indicates that wind speeds of 12.7 mph and greater,
blowing from the northwest, west and southwest occur 400 hours every year. Knowing the
number of hours the wind blows at different wind speeds the total cubic yards of sand
Table El:
Sand Transport Calculations
Existing
Yearly Australian Pine Stand Proposed Planting Plan
Wind Wind Speed Sand Accumulation Sand Accumulation
Speed Frequency Rate (yd3/yr) Rate (yd3/yr)
(mph) (hrs) VEG 1.5 VEG 2.0 VEG 1.5 VEG 2.0
12.7 144 64.8 64.8 64.8 64.8
13.8 108 64.8 64.8 10.8 10.8
15.0 60 51.0 51.0 12.0 12.0
16.1 36 39.6 39.6 9.0 9.0
17.3 28 40.6 40.6 8.4 8.4
18.4 8 14.4 14.4 3.2 3.2
19.6 12 27.0 21.0 5.4 5.4
24.2 4 19.6 15.2 3.8 3.8
S_321.8 311.4 117.4 117.4
trapped by the canopy can be estimated (see Table El). For this case it can be shown that
311.4 cubic yards of sand is trapped by the Australian pine stand each year. If these trees
were removed and the proposed planting plan installed, only 117.4 cubic yards of sand per
year would be trapped. Thus, the proposed plan is not adequate to maintain the existing
sand accumulation rate.
The VEG program can now be used to alter the proposed planting plan in order to
maintain the existing sand trapping capability. To do this, a trial and error method of
determining a more satisfactory proposed plant spacing and/or suggesting more
appropriate plant species capable of matching or exceeding the existing planting's sand
trapping rate is used. For this case, the VEG 1.5 program indicates that if the spacing of
the sand live oak was changed to 10 feet and planted at a height of 12 feet and the dune
prairie grass planted 1 foot apart at 1 foot height, then the planting plan would be
acceptable for its sand trapping capabilities.
Applying this case to VEG 2.0 yields additional information about the effects of the
vegetation on sand transport through and behind the plant canopy (see Figure E2). For
higher wind speeds, this program indicates a length of influence, L1, of 2,310 feet, a
deposition length, L2, within the vegetation of 300 feet, and an erosion length, L3, behind
the vegetation of 2,010 feet for the existing Australian pine stand. The proposed planting
plan gives an L1 of 946 feet, an L2 of 300 feet and an L3 of 646 feet. What this indicates is
that the Australian pine stand affects the wind structure and erosion behind the canopy for
a much longer distance than the proposed planting plan. With the asphalt drive directly
behind the vegetation, the actual length of sand erosion becomes 1,998 feet and 634 feet
behind the existing and proposed canopies, respectively (see Figure E3). The length of
erosion, L2, behind the canopy remains the same regardless of the ground surface material.
The only time that erosion actually occurs is when the surface is bare sand.
The above example demonstrates only one use of the VEG program. It is expected that it
will be used in numerous other applications including, the regulation of destructive coastal
dune activities and planning mitigation for destructive natural occurrences.
N
/ Developers Property
i .,i 
/ A
/ I
t00O
0 V
//
0 C
S300' 
/ .^C
Figure El. Plan View of Case Study Area.
LA
Deposition
Erosion
Figure E2. Crosssection showing lengths of Influence, Deposition and Erosion.
I O ,, 15, ,, / / Erosion
L L
L, L3.(Actual)
Figure E3. Crosssection of Example Problem showing lengths of Influence, Deposition
and Actual Erosion.
