Physical effects of vegetation on wind-blown sand in the coastal environments of Florida

Material Information

Physical effects of vegetation on wind-blown sand in the coastal environments of Florida final report
Sheppard, D. M ( Donald Max ), 1937-
Niedoroda, Alan W
University of Florida -- Coastal and Oceanographic Engineering Dept
Place of Publication:
Gainesville Fla
Dept. of Coastal and Oceanographic Engineering, University of Florida
Publication Date:
Physical Description:
1 v. (various pagings) : ill., charts ; 28 cm.


Subjects / Keywords:
Sand dunes -- Florida ( lcsh )
Beach erosion -- Computer programs -- Florida ( lcsh )
Revegetation -- Florida ( lcsh )
Coastal ecology -- Florida ( lcsh )
government publication (state, provincial, terriorial, dependent) ( marcgt )
bibliography ( marcgt )
non-fiction ( marcgt )


Includes bibliographical references.
General Note:
"February, 1992."
General Note:
"UF Project No. 4910451123312."
General Note:
"UPN No. 90100411."
General Note:
"Contract No. C-6851."
General Note:
This publication is being made available as part of the report series written by the faculty, staff, and students of the Coastal and Oceanographic Program of the Department of Civil and Coastal Engineering.
Statement of Responsibility:
principal investigator: D.M. Sheppard ; co-principal investigator: A.W. Niedoroda.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
002192803 ( ALEPH )
26521859 ( OCLC )
ALD2616 ( NOTIS )


This item has the following downloads:

Full Text

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 Wind-Blown 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

1. AGENCY USE ONLY (Leave bWank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVEREDFebruary 1992 Final Report
Physical Effects of Vegetation on Wind-Blown Sand C No. C-6851 in the Coastal Environments of Florida
6. A UT HOR(S)
D.M. Sheppard
A.W. Niedoroda
*University of FloridaGainesville, Florida 32611 PR 4910451123312
*Enviroanental Science & Engineering Inc. UPN No. 90100411 Gainesville, Florida
9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESSES) 10.SPONSORING /MONITORING U.S. Dept. of Commerce/NOAA Dept. of Eng. Reg. AGENCY REPORTNUMBER OCRM Coastal Management 1825 Connecticut A., N.W. 2600 Blair Stone Rd. Washington, D.C. 20235 Tallahassee, FL 32399
13. ABSTRACT (MaximumR200 woNUM)B
One and two-dimensinal, 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 coeffcients for crop plants similar in shape and size to the more common coastal vegetation in Florida.

Eolian sand transport / wind blown sand transport /
-canopy flow / sand trapping..

Standard Form 298 (Rev. 2-89) Prextrbed by ANSI Std. Z39-18 298_102

".. 7543-01-280-550

Partial funds for this project were provided by the Department of Environmental Regulation, 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.

Final Report

Physical Effects of Vegetation on Wind-Blown Sand
in the Coastal Environments of Florida
Principal Investigator: D.M. Sheppard
Department of Coastal and Oceanographic Engineering University of Florida
Co-Principal 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 vegetation 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 one-dimensional 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 two-dimensional 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 one-dimensional 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 program (Version 2.0) uses parameterized results from numerical experiments performed with the two-dimensional model thus some of the more important two-dimensional 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.

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 near-ground 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 one-dimensional 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 one-dimensional 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 two-dimensional model and program was developed using an extension of the turbulence algorithm used in the one-dimensional 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 one-dimensional model. A brief description of the capabilities and limitations of the one and two-dimensional 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 one-dimensional 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 two-dimensional model. The turbulence models used in both the one and two-dimensional 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 one-dimensional 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 one-dimensional 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 one-dimensional 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 two-dimensional 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 3090-600J). Details of the two-dimensional model, governing equations and the numerical algorithms used in the solution are given in Appendix C. The two-dimensional 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 two-dimensional 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 two-dimensional 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 beach-dune 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 one-dimensional 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 two-dimensional 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 two-dimensional 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 two-dimensional 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 two-dimensional model.
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
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, 81-99.
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." Boundary-Layer Meteorology, V.50, 227-275.
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." Boundary-Layer Meteorology, V.42, 293-311.
Baldocchi, D.D. and T.P. Meyers, (1988) "A Spectral and Lag-Correlation Analysis of
Turbulence in a Deciduous Forest Canopy." Boundary-Layer Meteorology, V.45,
Baldocchi, D.D. and T.P. Meyers, (1987) "Turbulence Structure in a Deciduous Forest."
Boundary-Layer Meteorology, V.43, 345-364.
Baldocchi, D.D. and B.A. Hutchison, (1987) "Turbulence in an Almond Orchard: Vertical
Variations in Turbulent Statistics." Boundary-Layer Meteorology, V.40, 127-146.
Baldocchi, D.D., S.B. Verma and N.J. Rosenberg, (1983) "Characteristics of Air Flow
Above and Within Soybean Canopies." Boundary-Layer Meteorology, V.25, 43-54.
Baldocchi, D.D., (1982) "Mass and Energy Exchanges of Soybeans: Microclimate-Plant
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, 1543-1544.
Eghlima, A. and C. Kleinstreuer, (1985) "Numerical Analysis of Attached Turbulent
Boundary Layers along Strongly Curved Surfaces." AIAA Journal, V.23, No.2,


Fast, J.D. and E.S. Take, (1988) "Evaluation of an Alternative Method for Numerically
Modeling Nonhydrostatic Flows over Irregular Terrain." Boundary-Layer Meteorology, 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, 285-304.
Finnigan, J.J. and P.J. Mulhearn, (1978) "Modelling Waving Crops in a Wind Tunnel."
Boundary-Layer Meteorology, V.14, 253- 277.
Finnigan, J.J., (1979) "Turbulence in Waving Wheat." Boundary Layer Meteorology,
V.16., 181-211.
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,
Garratt, J.R., (1990) "The Internal Boundary Layer A Review." Boundary-Layer Meteorology, V.50, 171-203.
Gong, W. and A. Ibbetson, (1989) "A Wind Tunnel Study of Turbulent Flow over Model
Hills." Boundary-Layer Meteorology, V.49, 113-148.
Halldin, S., (1985) "Leaf and Bark Area Distribution in a Pine Forest." The
Forest-Atmosphere Interaction, 39-58.
Haworth, D.C. and S.B. Pope, (1986) "A Generalized Langevin Model for Turbulent
Flows." Phys. Fluids, V.29, No.2, 387-405.
Hogstrom, U., (1988) "Non-Dimensional Wind and Temperature Profiles in the
Atmospheric Surface Layer: A Re-Evaluation." Boundary-Layer Meteorology, V.42,
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, 163-178.
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, 213-232.
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, 101-120.
Horikawa, K., S. Hotta and S. Kubota, (1982) "Experimental Study of Blown Sand on a
Wetted Sand Surface." Coastal Engineering in Japan, V.25, 177-195.


Hotta, S., N.C. Kraus and K. Horikawa, (1987) "Function of Sand Fences in Controlling
Wind-Blown Sand." Proc., Coastal Sediments '87, ASCE,, V.1, 12-14.
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, 1265-1281.
Hutchison, B.A., et. al., (1986) "The Architecture of a Deciduous Forest Canopy in
Eastern Tennessee, U.S.A." J. of Ecology, V.74, 635-646.
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, 1-88.
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,
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
Navier-Stokes Equations." J. of Computational Physics, V.59, 308-323.
Kimura, F. and R. Manins, (1988) "Blocking in Periodic Valleys." Boundary Layer Meteorology, V.44, 137-169.
Kleinstreuer, C., A. Eghlima and J.E. Flaherty, (1982) "New Higher-Order
Boundary-Layer Equations for Laminar and Turbulent Flow Past Axisymmetric
Bodies." AFOSR-TR-82-0588, 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, 729-741.
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, 611-628.
Li, Z., J.D. Lin and D.R. Miller, (1990) "Air Flow over and Through a Forest Edge: A
Steady-State Numerical Simulation." Boundary-Layer Meteorology, V.51, 179-197.
Maitani, T., (1989) "Wave-Like Wind Fluctuations Observed in the Stable Surface Layer
over a Plant Canopy." Boundary-Layer Meteorology, V.48, 19-31.


Mellor, G.L., (1973) "Analytic Prediction of the Properties of Stratified Planetary Surface
Layers." J. of the Atmospheric Sciences, V.30, 1061-1069.
Moeng, C.H., and J.C. Wyngaard, (1989) "Evaluation of Turbulent Transport and
Dissipation Closures in Second-Order 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 Second-Moment Closure." 27th Aerospace Sciences Meeting, January
9-12, 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, 548-553.
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,
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- ary-Layer Meteorology, V.43, 43-63.
Pinker, R.T., (1983) "The Canopy Coupling Index of a Tropical Forest." Boundary-Layer
Meteorology, V.26, 305-311.
Raynor, G.S., (1971) "Wind and Temperature Structure in a Coniferous Forest and a
Contiguous Field." Forest Science, V.17, No.3, 351-363.
Shaw, R.H. and I. Seginer, (1987) "Calculation of Velocity Skewness in Real and Artificial
Plant Canopies." Boundary-Layer Meteorology, V.39, No.4, 315-332.
Shih, T-H. and J.L. Lumley, (1986) "Second-Order Modeling of Near-Wall Turbulence."
Physical Fluids, V.29, No.4, 971-975.
Skupniewicz, C.E. R.F. Kamada and G.E. Schacher, (1989) "Turbulence Measurements
over Complex Terrain." Boundary-Layer Meteorology, V.48, 109-128.
Spalart, P.R., (1988) "Theoretical and Numerical Study of a Three-Dimensional
Turbulent Boundary Layer." J. of Fluid Mechanics, V.205, 319-340.
Speziale, C.G., (1989) "Discussion of Turbulence Modeling: Past and Future." Final
Report for Institute for Computer Applications in Science and Engineering, Report
No. 89-58, 23 p.

I .

Speziale, C.G., (1982) "On Turbulent Secondary Flows in Pipes of Noncircular
Cross-Section." International J. Engineering Sciences, V.20, No.7, 863-872.
Thom, A.S., (1971) "Momentum Absorption by Vegetation." Quart. J. R. Met. Soc.,
V.97, 414-428.
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
Plant-Air Layer." Bull. Nat. Inst. Agr. Sci., Ser. A, No.11, 19-65.
Wilson, J.D., (1988) "A Second-Order Closure Model for Flow Through Vegetation."
Boundary-Layer Meteorology, V.42, 371- 392.
Wilson, N.R. and R.H. Shaw, (1977) "A Higher-Order Closure Model for Canopy Flow."
J. Applied Meteorology, V.16, 1197-1205.
Wyngaard, J.C., (1987) "A Physical Mechanism for the Asymmetry in Top-Down and
Bottom-Up Diffusion." American Meteorological Society, V.44, No.7, 1083-1087.
Wyngaard, J.C. and R.A. Brost, (1984) "Top-Down and Bottom-Up Diffusion of a Scalar
in the Convective Boundary Layer." American Meteorological Society, V.41, No.1,
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 Post-Nourishment" 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, 529-538.


Appendix B
One-Dimensional Canopy Flow Model

Appendix B

One-Dimensional Canopy Flow Model
The current version of the FDNR VEG code employs a one-dimensional 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 one-dimensional 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 one-dimensional, second order closure model consists of five partial differential equations for the five variables i1, u'w', 'u', i' 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:
1) a+ & +CAU2 =,
at az
8u'w' 85 8 ( ___ q 8
2) +ww--2- qA + u'w' -Cq2- = 0,
tz z /z 3A2 Z
aii' au u) +2 q2 )+2q3 1
3) + 2i'w' k + ('u' -i -- +(-7=u, (
4 at 19Z a z 0Z 3A2 3 3As
_V__' ( -_'_' q ( 2 2
4) at az a z) 3A2 3'' +3 =3
)'w + -!'- (j q2 2q!3
5) 3 qA +2- w'w'- ) + =0,
at az 0Z 3A2 3 3As
where z is the height, t is time, A,, 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 = (7 + i7-V7 + j ) 1/2 .
The three length scales A,, A2 and A3 are defined in terms of a single length scale, i:
A, = 81t, A2 = 82f and A3 = 6,
where 81, 82 and S3 are constants. The length scale 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
a- = (2) ,Oz KZ
where u. is the friction velocity and x is von Karmon's constant. Also, the following relationships have been determined from measurements in a constant stress layer, Ulu' = -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:
82 = 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,
ffU1U 8-P-W
= 0, =0 and
= 0.
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,
8_11= 0,a 7 = 0
z z
= 0, and -v = 0
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).
Zn < U*c
= 22 (4)
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 one-dimensional 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:


a9 8F
-+- + 0 (5) t az
9 'u' a'u / I F -3qAj O9w'w'
_F 2u~
S 2w'u$ u +Cq
Equation 5 represents a system of five equations for five dependent variables contained in the vector 9. To obtain a numerical solution, the vertical domain is represented by discrete points, z;, which span from z0 to H. Between each grid point the spacing Az is defined as
Az+1/2 = zi1 z.
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 9 are then approximated by discrete values, 9, at each grid point i.
Using this approximate representation, the spatial derivative of the term F is approximated with finite differences:
~ ~ F (1+1/2) e(l-1/2) kZ z ~ Azj


F (+1) I (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, 9:
+F + 2 =0. (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.
____(n+i) (+) + F S(n+i) =0, (7) At
(n+i) n
and n indicates the time step. Equations 7 are nonlinear due to the terms ( and S (they require values for
Q+ which are yet unknown). Therefore, the terms F and S must be linearized in time:
(n+i) n F F + AQ and
a Q
(n+i) Sn
where and are the Jacobian matrices, herein denoted as G and H, respectively.
a Q a Q


Substituting the linearizations of F and into Equation 7, we get
A Q +'5b (GA Q, H Q =6JF + Sn
-At Z+H
which can be rewritten as
[I + At (6Gi + Hi)] A9 = At ( + ). (8)
This equation is a matrix system of equations for the dependent variables A 9 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 9 -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+1) 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 T! and the four Reynold's stresses u'w', v and iwl. The shear stress at the ground, which is required to determine the sand transport, can be obtained from the steady state solution as follows: = p2
where u. = V(-Ui'w)z1.
Here rg is the ground shear stress, p is the air density, u,. is evaluated at the ground (i.e. at zo). The ground shear stress r, is then used to determine the sand flux.



Part B

Model Parameters
The one-dimensional 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 61 for the plant types contained in the database are discussed below.
Leaf Area Density
The leaf area index (LAI) is a measure of the one-sided 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 non-coastal 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 B1-B10. 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 VI. The Cd and 81 for the reference plants were calibrated with the one-dimensional model and the measured wind velocity profiles obtained in the literature for the reference plants. The Cd and 81 were varied in the model to produce velocity profiles that best fit the existing measured wind velocity profiles. The parameter 81 was given just two values, 1.0 and 0.1, while the Cd was changed to produce the best fitting profile for each 81. The model profiles compared with the experimental profiles for each reference plant are shown in Figures B11-B20.


Table Bi:

NA = not available


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

Table B2:

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

1.20 H. = 26.3 ft.
L. = 26 ft.
1.00 0.80 =?0.60
0.40 0.20
0.00 o6""~~~o a
Leaf Area. Index

Figure B1

1.20 1.00 0.80 =!0.60
0.40 0.20 0.00 0.


,, Area ,,,.0I' Leaf Area Index

Figure B3

1.20 i H. = 3.9 ft.
L. = 2.5 ft.
1.00 0.80
0.40 0.20 S460.00~:o idb
Leaf Area Index
Figure B2

1.20 H. = 75.5 ft t. = 8.0 ft
1.00 0.80 20.60
0.40 0.20
0.0 .2 0.40" 6.0 ".0" 1.00
Leaf Area Index
Figure B4

H. = 4.6 ft L. = 1.3 ft

'''"... K6

Leaf Area Density Profiles of Reference Plants

1.00 0.40 0.20

1.20 H. = 13.9 ft
L. = 13.1 ft
0.00 . . . . . . . ,...............,......,
0.00 0.10 0.20 0... 6 .4 M 00,0.60
Leaf Area Index

Figure B5

1.20 H. = 2.5 ft
L. = 0.7 ft
1.00 0.80 =!0.60
0.00 2.00 4.00 6.00 8.00 10.00
Leaf Area Index

H. = 34.5 ft L. = 8.6 ft

. 1 r . .
0.00 0.20 0.4d 0.60 0.80 1.00 1.0
Leaf Area Index

Figure B6

1.20 H. = 3.5 ft
L. = 2.5 ft
*0 = 0.60
0.0 .
0.00 2. . :0%, 10.00
Leaf Area Index

Figure B7

Figure B8

Leaf Area Density Profiles of Reference Plants

1.20 H. =34.5 ft
L. =8.6 ft
0.00 '. .. "'. .
Leaf Area Index
Figure B9

1.20 H. = 4.1 ft L. = 0.6 ft
1.00 0.80 =0.60
0.40 0.20
0.00 0
Leaf Area Index
Figure B10

Experimental and Model Wind Velocity Profiles
of Reference Plants

+ //
+++++ experime d, =0.1 - d, = 1.0

....... ... ....... ....
2.00 4.00 6.00
Mean Velocity



2.00 1.50 1.00

2.00 1.50 CD
0.50 0.00

- ..
+++++ experimental 4 d, =0.1, C
4 1.0


0.20 0.40 0.60 0.50
Mean Velocity

1.00 1.,

Figure B11

Figure B12


+++++ experimental ---- d1=0.1.Cd
= 1.0, Cd

I 0.20


data = 0.15 = 0.17

2.00 1.50 .01.00

+ -
+++++ experime + d, = 0.1 d1 = 1.0



ntal data Cd = 0.05 C = 0.06

2.00 4.00 6.00 8.00 10.00
Mean Velocity

Figure B13

ntal data Cd = 0.18 ,Cd = 0.29
. . . ,


0.00 +,,

data = 0.06 = 0.065


0.50 0.00

0.40 0.60 0.80 1.00 1.20
Mean Velocity

Figure B 14

... a n.a........a...,.ac.c




and Model Wind Velocity Profiles of Reference Plants

Orng a

+ +++++ experim
d, =0.
- - d, = 1.



0.00 ec. i0t
Mean Velocity

Figure B15

u.00 1i1

+++++ expertmenfa
-- d, = 0.1, C --- d, = 1.0, C

I data d = 0.06 d = 0.075

1.00 2.00 3.00 4.0
Mean Velocity

Figure B16

+ /
+++++ experimental data
- d, =0.1, C, = 0.14 d, =1.0, Cd = 0.18
r. .....'... '' ' 2.h ' .00
0.50 1.00l ity .
Mean Velocity

1.00 0.50


/ / /

experimental data d1 = 0.1, Cd = 0.10 d1 = 1.0, Cd = 0.19

266 0 4.00 y
Mean Velocity

Figure B17

ental data .1, C, = 0.6 0, Cd = 0.5






Figure B 18

Experimental and Model Wind Velocity Profiles
of Reference Plants

+ -
+++++ experime -- d1= 0.1
d, =1.0

Mean Velocity



ntal data Cd = 0.5 Cd = 0.7
' .' 0

2.00. Wheat
+++++ experimental dat
- d = 0.1, C = d, = 1.0. Cd = 0.00 '' ' 66''',''6'' 0
Mean Velocity

Figure B20

2.00 1.50
~1.00 0.50-

20 CD



Figure B 19


Appendix C
Two-Dimensional Canopy Flow Model

Appendix C

Two-Dimensional Canopy Flow Model
A two-dimensional second order turbulence model for flow through a finite length canopy on a flat terrain has been developed in order to determine two-dimensional effects such as the dependence of the transition length on the canopy properties. Following is a description of the two-dimensional 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 two-dimensional second order turbulence model consists of six partial differential equations for the six dependent variables U, F, '7, u'u', v'v', and w'w'. They are:
a- + 0) = 0, (1) ax Oz
aU Ou2 o all OW- 6-U-, ap
+ + + + + = -CdA|ul, (2)
at x az 09x az ax
8-W a- r2 a-U-,& 8w'w' aP
-- + + + +z -CdA I, (3)
at ax + z
+ UW + w Cq2 + (-u~'' Cq2 + u'w'
_9 az ax~,auw' .k ~Oau'w a a aw'w9 ax ax az az) ax az 9Z ax
- iU/ i, (4)
__'_ au'u' aWU7 2) Oa
5t + U + W +2 W7 -Cq2 + 2'wat ax az Ox 9Z
3 a (qAj Ou'u' -a (qA7j' 2\ au'w' axk ax] a az ) az ax
= 3A 2 -3(3)


___ a 7i P7 (qA (qX a' (
at + ii + U7- ---q
ax az ax ax az az
q 27L and (6) 3A 3 3A3
aw'w' aw'w' _aw'w' a- 281
at1 w __az7 8W 2(W-7W-_C,2) az
+ 2 ax + 2u'w'-- '
- A. Sw'w' ..2 o1w'w' 8 2A'.
x ax) az ( az ) ax -az
3A2 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, 81, 82 and 83 are defined and discussed in Appendix A. These equations are simular to those of the one-dimensional turbulence model except that: 1) the equations for the Reynolds stress have advective and diffusion terms in two-dimensions, 2) there is an additional momentum equation for T7 (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 one-dimensional model with the addition of conditions for UT (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 two-dimensional second order turbulence model is a semi-implicit finite difference scheme. It is however, more complicated than the scheme employed in the one-dimensional model (Appendix A) due to two factors. First, since the equations are two-dimensional, a spatially factored scheme is necessary to efficiently solve the Reynolds stress and momentum equations. Secondly, since the momentum equations


contain the dynamic pressure term and the continuity equation is not automatically satisfied (as it is in the one-dimensional 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 two-dimensional 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 CI shows the grid system and variable location. Note that the Reynolds stresses R are defined at the nodes, the mean velocity components are "staggered" with the u;, and vij components defined on the vertical and horizontal cell faces respectively and the pressure Pi, is defined at the cell centers. Also defined on the nodes are values of the leaf area density Ai3. 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 4-7) for the stress components in vector form:
r-. 11- S2
+ -+ =-a-(8
lox Z(8)
7= (9)
uuT&W qAj1 -w-- Uu
uu'u'- 39u1
FV (10)
uv'W' qAj' 21


wu'w' 2qAi 8UZW- W'W1a7 GZ
wv q()
wwI'w 3qA 1
(wW-wI-UC2 2) 8 7 -78W
-2 Cq2(uu' - (12)
s 2 2 (iiw' Cq 2 Z)7 3 2 (3+3
2 =3 (12)
An implicit time marching algorithm (in delta form) is first introduced and can be represented as:
At'+ 56, (AARI) + s, (BA B) RA d = (RHS);1, (13) where AR3 = R,"+1 R," (n indicates time level), A and B are the Jacobian matricies
2 3
A = and B-=
and RHS2 is- G + .
13 13 13
The difference operators 6, g, 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.
R~ as
+ A (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 can then be obtained from
- Id+ (15) 13 i
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--ASA*-t(I +AtyB*) A R. =-t (~ 13
I+ Atb Aji At (IRt,;)AB =A (RHS)ij. (16)
Note that the difference between this equation (when expanded) and Equation 13 is the
a S
term At2(5xA -)8yB. This term is of order At2, thus only a small error is introduced
a R
by omitting it. Also, once the converged solution is obtained, AR3 = 0, RHSy = 0 and subsequently the error will also go to zero. Thus this error will not affect the steady-state solution.
Equation 16 can be written as a two step process: (I + At6xA) AR,AtRHS;1 and (17)
(I + AtvBi3) AR;, = ARj. (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 AR;3 (using current values for i and U).
3. Obtain Rn for A .
4. Update the boundary conditions.
5. Repeat steps 2-3 until AR3 -> 0.
Note that the term S,.. contains values of the mean velocity components. These values hi
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 1-3), 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:
-+ P-= -- -- -y;- CA|Iiiand
0 j ilw ww CdAlgiu.
-t+ Py --X CaAaiz Ot x Oz x Oz First by defining the following intermediate velocities W and U* as
-= -+P,and
at at Z'
and then applying standard finite difference approximations to Equations 19, 8 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:
+ O- = (2p 2P (21) t x Oz -t x Oz 5X2 + Oz2J
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 two-dimensional 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 two-dimensional model are shown in Figure C4. Upstream of the canopy a log-layer 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 one-dimensional 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 log-layer profile. In this region the bottom stress continually increases until the bare ground value is obtained.
Two-Dimensional Effects
As seen in the sample case above, a number of two-dimensional effects can occur which will affect the total transport of sand into and out of a canopy. The two-dimensional model has been used to predict the two-dimensinal 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 two-dimensional effects. The lengths are depicted in Figure C5 and are defined below:
L1: 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 one-dimensional model. However, if the canopy is shorter than the transition length then the minimum bottom stress will not decrease to the one-dimensional model value. The length L1 includes both the canopy length and the length of the wake region. Results from the two-dimensional model indicate that the wake region will range from 22 to 45 canopy heights, thus Li 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. L3-Lc, 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 two-dimensional model has been incorporated into a database for use in the FDNR VEG (Version 2.0) code. Results from the one-dimsensional 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 one-dimensional 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 two-dimensional database.


dx R+1-1
j-- xP dz
j+1j Figure C1 Staggard grid system and variable locations.
4.00 j



20 40 60

Figure C2 Canopy dimensions and coordinate system.




I .

2.0 1.5
1.0 0.5
0.0 -r,7
0.00 0.05 0.10 0.15 LAD Figure C3 Canopy leaf area density profile (LAD).

47 3

10 15
(at x/H = -20)



5 10 15 Velocity (at x/H = 0.5)


N 2-

Top of Canopy
--------------- ---------------


5 10 15 20 Velocity (at x/H = 30)

Top of Canopy
---------------- ------------




10 15 20 (at x/H = 45)


Figure C4 Velocity profiles for 2-D flow through canopy.

Top of Canopy
----------------- -----------

Top of Canopy
---------------------- --------




N 21-


n i


L1 deposition

Lc deposition
LL depositionAm L2

U 3


Figure C5 Definition of lengths L1,L2,L3 and Lc.

erosion -

erosion \I

Appendix D
Flow Over an Irregular Boundary

Appendix D

Flow Over an Irregular Boundary
The first step in the development of a two-dimensional 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 two-dimensional flow through canopies on flat terrain. The flat terrain problem was addressed using the two-dimensional model developed in Appendix C. The approach taken here to develop a two-dimensional 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 C1-C7) 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 two-dimensional canopy flow model for genereal 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 Navier-Stokes equations:
-+-x 0 (1) Ox Oz'
S--+ -+ + = V + a2U) and (2) t Ox Oz OZ2
GQ -W m ft7&W2 ap a2 2U7 + '+ j+ -- V + (3) 8t a z oZ Wx2 aZ2These equations are similar to equations C1-C3 but with the Reynolds stress terms replaced with viscous stress terms. In developing a solution procedure, these equations



have been transformed onto a general curvaliner coordinate system and solved using finite volume methods. The numerical algorithm is essentially similar to the time-split 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 efficency 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 shore-normal by 1500 feet shore-parallel. 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 non-native 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 one-dimensional 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 two-dimensional 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

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
L _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.



/ / / /
/ Developers Property
/ I / I
/ V
/ ,?~ ,Ij i?;
/ ~e C -~ C~
/ E~ 'P
// C. ~ C/ 4'~ 'P 'P
0 / (fl
/ ~0
/ ~y, ~
- ,~ C.
'P 'P 'P 2~.
o /
o / *~, ',, q~,
// .~
/ / /
/ 'P 'P 'P
/ /
/ Vg~ V~ 'P
/ c~
/ 'P 'P 'P
C. C. ~ ~p
Vg~ 'P
// I <-30

Figure El. Plan View of Case Study Area.


Deposition Erosion
F/ Figure E2. Cross-section showing lengths of Influence, Deposition and Erosion.





Figure E3. Cross-section of Example Problem showing lengths of Influence, Deposition
and Actual Erosion.