Ecological Archives E085-101-A2

Merel B. Soons, Gerrit W. Heil, Ran Nathan, and Gabriel G. Katul. 2004. Determinants of long-distance seed dispersal by wind in grasslands. Ecology 85:3056–3068.

Appendix B. The Markov chain Synthetic Turbulence Generation (STG) model. A pdf file of this appendix is also available.

The Markov chain STG model is a version of the forest seed dispersal model of Nathan et al. (2002b), adapted for grassland ecosystems. The general structure is presented in the Methods section and by Katul and Albertson (1998) and Katul and Chang (1999). In this appendix the estimation of the Eulerian flow statistics needed to drive the Lagrangian dispersion is described, followed by descriptions of the computational grid, numerical scheme, boundary conditions and closure constants.

Upon temporal and spatial averaging the conservation of momentum equations, and following the closure approximations of canopy flows, the Wilson and Shaw (1977) closure model simplifies to the following set of ordinary differential equations (ODEs):

 Mean Momentum: 
Tangential Stress Budget: 
Longitudinal Velocity Variance:
Lateral Velocity Variance: 
Vertical Velocity Variance:


Where is a characteristic velocity scale, Cd is a drag coefficient, a(z) is a leaf area density, j = aj × Lws , with Lws a characteristic length scale specified using the formulation given by Katul and Albertson (1998) and not permitted to increase at a rate larger than k, and a1, a2, a3 and C are determined so that the flow conditions well above the canopy reproduce well-established surface layer similarity relations. With estimates of the five constants (a1, a2, a3, C, and ), the above five ODEs can be solved iteratively for the five flow variables , , , , and , which are used to drive the Lagrangian model.


The computational grid

The computational flow domain was set from zero to 20 × h. The grid node spacing is z = 0.005 m. This grid density was necessary due to rapid variability in leaf area density close to the canopy top. Parameter values at the exact location of the seed are calculated by interpolation between the grid nodes, or extrapolation in the case of x3 > 20 × h. To ensure seeds do not exit the atmospheric boundary layer during the computation of dispersal trajectories, the vertical position of a seed is not allowed to increase above 800 m.

The numerical scheme

The five ODEs for the Wilson and Shaw (1977) model were first discretized by central differencing all derivatives. An implicit numerical scheme was constructed for each ODE with boundary conditions to be discussed in the following section. The tridiagonal system resulting from the implicit forms of each discretized equation was solved using the Tridag routine from Press et al. (1992; pp.42–43) to produce the turbulent statistic profile.  Profiles for all variables were initially assumed, and a variant of the relaxation scheme described by Wilson (1988) was used for all computed variables. Relaxation factors as small as 5% were necessary in the iterative scheme because of the irregularity in the leaf area density profile. The measured leaf area density was interpolated at the computational grid nodes by a cubic-spline discussed in Press et al. (1992; pp.107–111) to insure finite second derivatives of a(z). Convergence is achieved when the maximum difference between two successive iterations in  did not exceed 0.0001%. We checked that all solutions were independent of z (as described in Katul and Albertson 1998). Calculation of dt is described in Appendix A.

Boundary conditions and closure constants

Typically, the well-established flow statistics in the atmospheric surface layer provide convenient upper-boundary conditions for closure models.  The boundary conditions used are:


Where is the standard deviation of any flow variable  (= ), Au = 2.2, Av  = 2.0, and Aw = 1.4 (Panofsky and Dutton 1984).

The closure constants are dependent on the choice of the boundary conditions and are determined by assuming that in the atmospheric surface layer (z > 2h), the flux-transport term is negligible and that , , and become independent of z for near-neutral conditions. These simplifications result, after some algebraic manipulations (e.g., see Katul and Albertson 1998; Katul and Chang 1999), in the following relationships between Au, Av, and Aw and a2, a3, and C:


Where Aq = (Au2 + Av2 + Aw2)1/2. The closure constant a1 is determined by noting that the eddy-diffusivity is k × (z – d) × u* in the surface layer. Hence, q1 becomes identical to k × (z – d) × u* leading to a1 = 1/Aq. The above equations are the first analytic expressions relating closure constants to ASL boundary conditions for the Wilson and Shaw (1977) model as described by Katul and Albertson (1998) and Katul and Chang (1999). Table B1 summarizes the closure constants used resulting from our choice of Au, Av, and Aw.


Table B1. Closure constants used in the Markov chain STG model for Au  = 2.2, Av = 2.0, and Aw = 1.4.

Closure constant
(Wilson and Shaw 1977)















Katul, G. G., and J. D. Albertson. 1998. An investigation of higher order closure models for a forested canopy. Boundary-Layer Meteorology 89:47–74.

Katul, G.G., and W.H. Chang. 1999. Principal length scales in second-order closure models for canopy turbulence. Journal of Applied Meteorology 38:1631–1643.

Nathan, R., G.G. Katul, H.S. Horn, S.M. Thomas, R. Oren, R. Avissar, S.W. Pacala, and S.A. Levin. 2002b. Mechanisms of long-distance dispersal of seeds by wind. Nature 418:409–413.

Panofsky, H. and J. Dutton. 1984. Atmospheric Turbulence: Models and Methods for Engineering Applications. John Wiley, New York.

Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery. 1992. Numerical Recipes in Fortran. Cambridge University Press.

Wilson, J.D. 1988. A second order closure model for flow through vegetation. Boundary-Layer Meteorology 42:371–392.

Wilson, N.R., and R.H. Shaw. 1977. A higher order closure model for canopy flow. Journal of Applied Meteorology 16:1198–1205.

[Back to E085-101]