## Abstract

Triaxial sonic anemometer velocity measurements vertically arrayed at six levels within and above a pine forest were used to examine the performance of two second-order closure models put forth by Wilson and Shaw and by Wilson. Based on these measurements, it was demonstrated that Wilson’s model reproduced the longitudinal velocity standard deviation *σ*_{u} better than did Wilson and Shaw’s model. However, Wilson and Shaw’s model reproduced the measured mean velocity 〈*u*〉 near the forest–atmosphere interface better than Wilson’s model did. The primary mechanisms responsible for discrepancies between modeled and measured 〈*u*〉 and *σ*_{u} profiles were investigated.

The conceptual formulations of these two closure models differ in the characteristic length scales and timescales used in the closure parameterizations of the mean turbulent kinetic energy dissipation rate term, the pressure–strain rate term, and the flux-transport term. These characteristic length scales were computed and compared with measured integral length scales inside the canopy. A discussion on how these length scales compare with the mixing layer analogy also is presented.

## Introduction

The quest for simple yet robust models of turbulent velocity statistics is motivated, in part, by the need to infer scalar sources, sinks, and turbulent fluxes at multiple levels within the canopy sublayer (CSL) of vegetated systems (e.g., Lewellen et al. 1980; Meyers and Paw U 1986, 1987; Leclerc et al. 1988; Raupach et al. 1992; Meyers and Baldocchi 1991; Baldocchi 1992; Katul et al. 1997b). Second-order closure models provide a practical framework for computing needed velocity statistics for modeling scalar transport within the CSL. Such models alleviate some of the theoretical objections to first-order closure and other phenomenological models (Albini 1981; Raupach and Thom 1981; Li et al. 1985; Lee et al. 1994; Wilson et al. 1998; Wilson 1989). In addition, their ability to reproduce measured second-order moments is comparable to third-order closure models for the CSL (see, e.g., Katul and Albertson 1998). As noted by Launder (1996), greater model complexity does not necessarily enhance model performance for many turbulent flow types, including CSL flows (Raupach 1988). Thus far, two second-order closure models have been proposed and tested successfully in the CSL of short vegetation (Wilson and Shaw 1977; Wilson 1988). However, their relative performance for tall forests with complex leaf area density distribution was not rigorously assessed.

In CSL second-order closure modeling, the conservation of streamwise mean momentum coupled with budget equations for tangential stresses results in three terms that require closure modeling, namely, the pressure–strain term, the momentum flux-transport term, and the turbulent kinetic energy dissipation rate term. Modeling these three terms requires specification of an empirical characteristic length scale or timescale as a function of the flow variables being modeled or computed (hereinafter referred to as principal scale). In the Wilson and Shaw (1977) closure model, an empirical length scale that varies linearly with height within the canopy (unless restricted by foliage drag and local leaf area density) was proposed based on arguments and measurements by Seginer et al. (1976). Alternatively, the Wilson (1988) model, which is based on Launder et al.’s (1975) and Hanjalic et al.’s (1980) approach, proposes a relaxation timescale defined by the ratio of turbulent kinetic energy *K* and mean kinetic energy dissipation rate ɛ. The arguments and justification for these scales are empirical; yet, they define principal length scales or timescales for important redistribution, transport, and dissipation mechanisms within the CSL. At present, there is no basis to judge the generality of either formulation for in-canopy turbulence. Differences in the turbulence scales translate to differences in model performance between the two closure schemes and have motivated a comparison between Wilson and Shaw’s (1977) and Wilson’s (1988) closure models for a tall uniform pine forest with a complex leaf area density profile.

Specifically, the objective of this study is to compare field measurements with modeled mean velocity, turbulent stress, and longitudinal and vertical velocity variances by these two closure models. Attention is devoted to similarity (or dissimilarity) in principal length scale (or timescale) variation within the canopy for the two closure models. How these length scales relate to the mixing layer analogy that was proposed recently by Raupach et al. (1996) also is discussed. The outcome of this comparison will guide modeling of the velocity statistics within the canopy for use in scalar transport models (Raupach 1989; Katul et al. 1997b). Furthermore, the ability of such closure models to predict third-order moments that potentially may be used in random-flight particle trajectory models also is considered.

## Theory

For a steady-state, neutral flow, the time- and horizontally averaged equations for momentum and Reynolds stresses are

where *x*_{i} (*x*_{1} = *x, x*_{2} = *y, x*_{3} = *z*) are the longitudinal, lateral, and vertical directions, respectively; *u*_{i} (*u*_{1} = *u, u*_{2} = *υ, u*_{3} = *w*) are the instantaneous velocity components along *x*_{i}; *p* is the static pressure normalized by the air density; *ν* is the kinematic viscosity; the overbar and 〈 〉 denote time and horizontal averaging, respectively; and primes and double primes denote departures from the temporal and horizontal averaging operators, respectively, as discussed in Raupach and Shaw (1982). All double-primed terms in (1) arise because of the horizontal averaging of the multiply connected air spaces within the CSL; hence, they are directly related to canopy-foliage effects on the flow statistics.

### Second-order closure approximations of Wilson and Shaw (1977)

In Wilson and Shaw (1977), the closure approximations for the streamwise momentum equation are

Such an approximation assumes that the form drag by the canopy is much larger than the viscous drag, and it employs *C*_{d} as a drag coefficient and *a*(*z*) as the plant area density at height *z.* The closure for the Reynolds stress equations are similar to those proposed by Mellor (1973):

where *q* = 〈*u*^{′}_{i}*u*^{′}_{i}〉^{1/2} is a characteristic turbulent velocity; 〈ɛ〉 is the mean rate of viscous dissipation; *λ*_{1}, *λ*_{2}, and *λ*_{3} are characteristic length scales for the triple-velocity correlation, the pressure–velocity gradient correlations, and the viscous dissipation, respectively; *δ _{ik}* is the Kronecker delta function; and

*C*is a similarity constant to be determined. As with typical second-order closure models, triple-velocity products are closed using a gradient–diffusion approximation. Pressure–velocity gradients are modeled on return-to-isotropy principles following Mellor (1973), Mellor and Yamada (1974), and Donaldson (1973); and viscous dissipation is assumed to be isotropic and dependent on the local turbulence intensity (Mellor 1973). Upon replacing these approximations in (1), and assuming horizontal homogeneity, the Wilson and Shaw (1977) model simplifies to

with

where *k*_{υ} (=0.4) is von Kármán’s constant; *α* is a constant to be determined experimentally; and *a*_{1}, *a*_{2}, *a*_{3}, and *C* are determined such that the flow conditions well above the canopy reproduce established surface layer similarity relations (e.g., Monin and Yaglom 1971). The term *L*_{WS}(*z*) is the empirical principal length scale of the Wilson and Shaw (1977) closure model (see also Lee et al. 1994) and is not permitted to increase at a rate larger than *k*_{υ} with increasing depth. In (4f), *d* is the zero-plane displacement height determined from the“center of pressure” method described in Thom (1971), Shaw and Pereira (1982), and Paw U and Meyers (1989) and is given by

where *h* is the canopy height. In the original formulation of Wilson and Shaw (1977), the zero-plane displacement height correction was not accounted for explicitly. The zero-plane displacement height correction is necessary because, in the atmospheric surface layer (ASL), the mixing length above vegetated surfaces is *k*_{υ}(*z* − *d*).

With estimates of the five constants (*a*_{1}, *a*_{2}, *a*_{3}, *C,* and *α*), the five ordinary differential equations in (4a)–(4e) can be solved iteratively for the five flow variables 〈*u*〉, 〈*u*′*w*′〉, 〈*u*′^{2}〉, 〈*υ*′^{2}〉, and 〈*w*′^{2}〉. The zero-plane displacement depends on 〈*u*′*w*′〉 and must be included in the iterative solution as an additional unknown. Details of the computational technique, determination of closure constants, and boundary conditions are presented in the appendix. The closure constants used in the Wilson and Shaw (1977) model are presented in Table 1. For notational simplicity, hereinafter the operator 〈 〉 defines simultaneous temporal and horizontal averaging.

### Wilson’s (1988) closure model

The Wilson model distinguishes between low-frequency-band shear production and high-frequency-band wake production of turbulence. The timescale for the low-frequency band is the ratio of mechanical production *K* and the dissipation rate ɛ. For such a band, the budget equations for the covariance tensor of the low-frequency band (surface kinetic energy band) are written as

The *P*_{ij} are the mechanical production terms and are given by

The *R*_{ij} are the redistribution terms, given by

(*ρ*_{0} is the air density) and parameterized as

The components of (7) are defined as

where

*n*_{k} = *δ*_{k3} is the unit-normal vector (in the vertical direction); *n _{i}*,

*n*, and

_{j}*n*are unit-normal vectors;

_{m}*δ*is the Kronecker delta function;

_{ij}*c*

_{1},

*c*

_{2}, and

*c*

^{′}

_{1}are similarity constants; the function in the bracketed right-hand side of

*β*

_{1}is the wall function, whose value is unity in the ASL (Wilson 1988; Gibson and Launder 1978; Shir 1973);

*K*is the mean turbulent kinetic energy; ɛ is the mean mechanical kinetic energy dissipation rate (to be formulated below); and

*A*

_{q}=

*q*/

*u*∗, where

*q*is determined from ASL similarity relations, and

*u*∗ is the friction velocity above the canopy. In traditional second-order closure models, (10a) and (10b) are added to account for the fact that a rigid wall modifies the pressure field, thus impeding the transfer of energy from the streamwise direction to that direction normal to the wall. The relative magnitude of the shear stress also is diminished by this term (Rosten and Worrell 1988; Shir 1973). Wilson (1988) adopted this form without any modification for CSL flows. Equations (8) and (9) represent the return-to-isotropy and fast-response contributions to the Reynolds stress budget. The fast-response term is associated also with coherent motion for CSL flows as demonstrated by Zhuang and Amiro (1994).

The *T*_{ij} in (5) are the turbulent transport terms, de- fined as

This term is modeled after Hanjalic and Launder (1972) as modified by Wilson (1988) for canopy flows via the constant *c*_{s} and is given by

where

LAI is the total leaf area index; *h* is, as before, the canopy height; and *c*_{so} is a similarity constant (see Table 1) whose value is derived from measurements by Hanjalic and Launder (1972).

The ɛ_{ij} terms in (5) are the low-frequency shear kinetic energy dissipation rate tensor and are given by

with

where ɛ_{cc} is the correct form of the mean kinetic energy dissipation rate in the absence of a canopy. In all of the above equations, *z* is replaced by (*z* − *d*) when *z* is greater than *h* to account for the zero-plane displacement height above vegetation. The zero-plane displacement height was not considered in the original model of Wilson (1988).

Upon replacing (6)–(13) in (5), the Wilson (1988) model, after some algebraic manipulations, simplifies to the following equations:

where

With estimates of the five constants (*c*_{1}, *c*_{2}, *c*^{′}_{1}, *c*_{3}, and *c*_{so}), the five ordinary differential equations (14a)–(14e) are solved iteratively for the flow variables 〈*u*〉, 〈*u*′*w*′〉, 〈*u*′^{2}〉, 〈*υ*′^{2}〉, and 〈*w*′^{2}〉 if appropriate boundary conditions are specified. These constants commonly are determined by matching the closure model to well-known ASL flow properties (see appendix for further discussion). As in the Wilson and Shaw (1977) model, the zero-plane displacement height, defined in (4g), must be included in the iterative solution for the Wilson (1988) model.

### Boundary conditions

Typically, the flow statistics in the ASL provide convenient upper boundary conditions for CSL closure models. For both models, the boundary conditions used are

where *σ*_{θ} is the standard deviation of any flow variable *θ* (=〈*θ*′^{2}〉^{1/2}); *A*_{u} = 2.4, *A*_{υ} = 1.95, and *A*_{w} = 1.25 are derived from many ASL field experiments for near-neutral atmospheric stability conditions (see Panofsky and Dutton 1984, p. 160); and *A*_{m} (=5.5) at *z*/*h* = l.44 is determined from sonic anemometer measurements as described in the next section. Table 1 summarizes the closure constants used by both models resulting from our choice of *A*_{u}, *A*_{υ}, and *A*_{w}.

## Experimental setup

Much of the experimental setup is described in Katul and Albertson (1998); however, for completeness, a brief summary is presented. The experiment was carried out from 25 May until 11 June 1997 at the Blackwood division of the Duke Forest near Durham, North Carolina (36°2′N, 79°8′W; elevation = 163 m above mean sea level). The site is an even-aged southern loblolly pine stand, established from seedlings planted at 2.4 m × 2.4 m spacing in 1983 following clear-cutting and burning (Ellsworth et al. 1995; Katul et al. 1997a,b). The canopy height neighboring the 20-m tall aluminum walk-up tower is 14 m (±0.5 m). The three velocity components *u*_{i} and virtual potential temperature *T*_{a} were measured at six levels using five Campbell Scientific, Inc., CSAT3 triaxial sonic anemometers and a Gill Instruments, Ltd., Solent sonic anemometer. The CSAT3 anemometers were positioned at *z* = 4.15, 5.95, 9.65, 13.16, and 15.91 m above the ground surface. One Solent anemometer was mounted on a short boom attached to the tower, at 20.70 m above the forest floor surface. A detailed comparison between the CSAT3 and the Solent sonic anemometers was performed on 11–13 October 1997 at *z*/*h* = 1.14. No systematic biases in measured 〈*u*〉, 〈*u*′*w*′〉, 〈*u*′^{2}〉, 〈*υ*′^{2}〉, and 〈*w*′^{2}〉 between the two sonic designs were found (see Katul and Albertson 1998). All signals from the six sonic anemometers were sampled at 10 Hz using a Campbell Scientific, Inc., CR9000 datalogger. The experiment produced 171 runs, with each run representing a 30-min sampling duration (18 000 measurements per flow variable per run). Of the 171 runs, 64 runs were collected during nighttime conditions and were excluded from the current analysis. The experiment resulted in 69 runs that satisfied three criteria:

the mean wind direction was between 110° and 250° (minimizing tower interference for the CSAT3 anemometers located just above the canopy),

the standard deviation of the wind direction did not exceed 50°, and

the

*u*∗ measured by the Solent anemometer exceeded 0.15 m s^{−1}.

The selected 69 runs form an ensemble at each height to which closure model predictions will be compared. To determine *A*_{m} at *z*/*h* = 1.44 from the 69 runs, a regression analysis (〈*u*〉 = *A*_{m}*u*∗) using the Solent measurements was performed, resulting in *A*_{m} = 5.5. The shoot silhouette area index, a value analogous to the LAI, was measured in the vertical at about 1-m intervals by a pair of LI-COR LAI-2000 plant canopy analyzers on 4 June 1997. The overall leaf area for this stand is 3.8 m^{2} m^{−2} (see Table 2 for detailed measurements of leaf area density).

## Results and discussion

In this section, comparisons between the two second-order closure model predictions and measurements are presented. A discussion on differences in characteristic length scales for the two closure models also is presented.

### Comparisons with measurements

#### First and second moments

In Fig. 1, the ensemble averages of the measured profiles of normalized *a*(*z*), 〈*u*〉, 〈*u*′*w*′〉, 〈*u*′^{2}〉^{1/2}, 〈*w*′^{2}〉^{1/2}, and *K* are shown with bounds of one standard deviation. As in all canopy sublayer field experiments, it is assumed that the collection of time averages, when properly normalized, forms an ensemble of realizations. For convenience, the numerical values of the measured ensemble mean and one standard deviation for relevant velocity statistics are presented in Table 3. The modeled profiles from the two closure models also are shown and compared to the measurements in Table 4. Table 4 reports the root-mean-square error (rmse) and regression statistics between modeled and measured flow variables. From Fig. 1 and Table 4, the Wilson and Shaw (1977) model (dashed line) reproduced the measured 〈*u*〉 better than did Wilson’s (1988) model (solid line); however, the Wilson (1988) model reproduced the *σ*_{u} profile better than the Wilson and Shaw (1977) model did. The Wilson (1988) model reproduced the measured 〈*u*′*w*′〉 profile better than did the Wilson and Shaw (1977) model in the lower layers of the canopy. All other modeled flow variables agree well with the measurements.

In addition, the computed zero-plane displacement heights by the two models (see Table 4) agreed with the widely used ⅔ value for natural surfaces (Brutsaert 1982, p. 116). Using the computed *d* and the measured *A*_{m}, an estimate of the momentum roughness height (*z*_{0}) can be obtained for this forest stand following algebraic manipulation of the logarithmic velocity profile formulation to give *z*_{0} = (*z* − *d*) exp(−*k*_{υ}*A*_{m}). With *z* = 1.44*h, d* ≈ (⅔)*h,* and *A*_{m} = 5.5, the resulting *z*_{0} = 0.09*h,* which is in agreement with the widely used *z*_{0} = 0.1*h* relationship (Brutsaert 1982).

Given the uncertainty in lower boundary conditions, a sensitivity analysis was performed in which gradients of zero for 〈*u*′*w*′〉, 〈*u*′^{2}〉^{1/2}, and 〈*w*′^{2}〉^{1/2} were specified at *z*/*h* = 0. It was found that at *z*/*h* > 0.1 the influence of the lower boundary conditions was dissipated and the two solutions (i.e., zero flux and zero state) converged.

To understand the mechanisms responsible for discrepancies between model calculations and measurements of 〈*u*〉 (by the Wilson model) and *σ*_{u} (by the Wilson and Shaw model), two analyses were performed. The first considers the third-moment closure models, while the second additionally considers components of the turbulent kinetic energy budget.

#### Third moments and closure approximations

In Fig. 2, ensemble averages of measured 〈*w*′*u*′*u*′〉, 〈*w*′*w*′*w*′〉, and 〈*w*′*u*′*w*′〉 profiles are shown with bounds of one standard deviation. It is evident that the predictions from both models compared poorly with measurements (see Table 4 for quantitative comparisons). What is alarming from the comparison in Fig. 2 is that the Wilson and Shaw (1977) model incorrectly predicted both magnitude and sign of 〈*w*′*u*′*u*′〉 close to the canopy–atmosphere interface (*z*/*h* = 0.95–1.3). This term is central to proper modeling of *σ*_{u}.

To demonstrate explicitly the role of 〈*w*′*u*′*u*′〉 in modeling *σ*_{u}, the Wilson and Shaw (1977) budget for 〈*u*′*u*′〉 can be rearranged so that

Since *K* (and hence *q*), 〈*u*〉, and 〈*u*′*w*′〉 profiles compare favorably with measurements, the only “candidate” term producing *σ*_{u} overestimation in Fig. 1 is *d*〈*w*′*u*′*u*′〉/*dz.* To test this hypothesis, consider the following in Fig. 2.

The measured 〈

*w*′*u*′*u*′〉 does not change much with increasing*z*for*z*/*h*= 0.6–0.9 (in fact, the measurements show a mild increase in 〈*w*′*u*′*u*′〉 with increasing*z*).In contrast, modeled 〈

*w*′*u*′*u*′〉 rapidly decreases with increasing*z*for*z*/*h*= 0.6–0.9.

This unrealistically large and negative *d*〈*w*′*u*′*u*′〉/*dz* leads to an unrealistic increase in 〈*u*′*u*′〉 as evidenced by (16) and the results in Fig. 1. Hence, the failure of the Wilson and Shaw (1977) model to reproduce measured *σ*_{u} for *z*/*h* = 0.6 and 0.9 can be attributed to poor parameterization of *d*〈*w*′*u*′*u*′〉/*dz.* In fact, the computed *σ*_{u} by the Wilson and Shaw (1977) model for their corn canopy (see their Fig. 3) also overestimated their measured *σ*_{u} for the *z*/*h* = 0.6–1.0 layers.

All in all, while the agreement between model predictions and measurements is poor, the scatter in the measured third moments is large (at least when compared to the scatter in measured second moments). As discussed by Shen and Leclerc (1997), thermal stratification impacts third moments much more than second moments, implying that *u*∗ and *h* are not the only normalizing variables for third moments. Hence, in addition to the gradient–diffusion theory failure to describe third moments from second moments, the neutral flow approximation further compounds this error in the closure schemes that parameterize third moments from second moments.

#### Turbulent kinetic energy dissipation rates

Having demonstrated why the Wilson and Shaw (1977) model overestimates 〈*u*′*u*′〉, we investigate next why the Wilson (1988) model overestimates 〈*u*〉 between *z*/*h* = 0.75–1.3. Rearrangement of (14b), after some algebraic manipulations, produces

Since *K,* 〈*w*′*w*′〉, and 〈*w*′*u*′〉 are reproduced well by the Wilson (1988) model (see Fig. 1), the two remaining terms that may explain overestimation in 〈*u*〉 are ɛ and *d*〈*w*′*u*′*w*′〉/*dz.* To examine which of these two parameters is more critical, two regions of the flow domain are considered: *z*/*h* = 0.75–1.0 and *z*/*h* > 1.0. In the region *z*/*h* = 0.75–1.0, it is clear that modeled *d*〈*u*〉/*dz* is larger than measured *d*〈*u*〉/*dz* (leading to an overestimation in 〈*u*〉). In the region *z*/*h* > 1.0, modeled *d*〈*u*〉/*dz* is much smaller than measured *d*〈*u*〉/*dz.* Using (17) and Figs. 1 and 2 one can consider the following points.

In the region

*z*/*h*= 0.75–1.0 (see Fig. 2), the measured*d*〈*w*′*u*′*w*′〉/*dz*is negligible. The*d*〈*w*′*u*′*w*′〉/*dz*predicted by the Wilson (1988) model is large in magnitude and negative in sign. Such a computed 〈*w*′*u*′*w*′〉 gradient leads to a large positive*d*〈*u*〉/*dz,*as evidenced by (17), with all else being modeled well (recall that*c*_{2}= 0.6 < 1). Hence, from this argument, it is clear that the Wilson (1988) modeled 〈*u*〉 is sensitive to the 〈*w*′*u*′*w*′〉 parameterization.In the region

*z*/*h*= 1.0–1.3 (see Fig. 2), both measured and modeled*d*〈*w*′*u*′*w*′〉/*dz*are small. Hence, the discrepancy between measured and modeled*d*〈*u*〉/*dz*cannot be attributed to the 〈*w*′*u*′*w*′〉 closure model (as before). The only remaining variable in (17) is ɛ. A negligible*d*〈*w*′*u*′*w*′〉/*dz*coupled with a small ɛ leads to a small*d*〈*u*〉/*dz*by (17) (in agreement with the comparisons of Fig. 1). Whether the computed ɛ by the Wilson (1988) model underestimates the actual dissipation rate cannot be tested directly because no direct dissipation rate measurements were available. An indirect test of this underestimation would be a comparison between the dissipation rates computed by the Wilson and Shaw (1977) and the Wilson (1988) models. In Fig. 3, the dissipation rates calculated by both models are compared. Notice that the dissipation rates calculated by Wilson’s (1988) model are consistently smaller than those from Wilson and Shaw’s (1977) models, in qualitative agreement with the above hypothesis as to why*d*〈*u*〉/*dz*is small for*z*/*h*> 1.0. In addition, a sensitivity analysis was performed to ensure that the above finding is not due to our estimate of*C*_{d}. We varied*C*_{d}from 0.15 to 0.30 and found that computed 〈*u*〉 systematically overestimates measured 〈*u*〉 at*z*/*h*= 1.1 for all*C*_{d}∈ [0.15–0.30].

### Principal length scale comparisons

As stated earlier, there is no basis to judge the generality of either the length scale or timescale formulations for in-canopy turbulence closure models. As evidenced from (4a)–(4e) and (14a)–(14e) these length scales and timescales influence important redistribution, transport, and dissipation mechanisms within the CSL. Given the good model agreements with measurements, it is important to determine how different these length scales and timescales are inside the canopy. From Wilson’s (1988) model, the low-frequency timescale is governed by the available turbulent kinetic energy (=*K*) and the mean dissipation rate (=ɛ). To compare the Wilson (1988) timescale with Wilson and Shaw’s (1977) empirical length scale, these two variables can be combined via dimensional analysis to define a principal length scale (*L*_{W88}) using

where *γ*_{W} is a constant. In order for the *L*_{W88} length scale to match its expected value in the neutral ASL [i.e., *k*_{υ}(*z* − *d*)], this constant must be

Because of the availability of simultaneous 10-Hz time series measurements for all 69 runs at six levels, both calculated length scales (*L*_{W88} and *L*_{WS}) can be compared to measurements within and above the canopy. The measured length scale was determined from *K* and the integral timescale *I*_{3} using

where *I*_{3} is determined from the autocorrelation function using

where *s* is time lag and *t* is time. For each of the 69 runs, the integration in (21) is terminated after the first zero crossing as discussed in Katul et al. (1997b). Since the integral timescale *I*_{3} measures a decorrelation time and hence is related to the dissipation time (or loss of coherence) of a characteristic eddy with mean kinetic energy *K,* then *I*_{3} is an indicator of a relaxation time. For consistency with Wilson’s (1988) low-frequency timescale argument, *K*^{1/2} (vis-à-vis 〈*u*〉) was chosen as the characteristic velocity to convert *I*_{3} to length as in (20). Note that Wang and Takle (1995) demonstrated that *K*^{1/2} and 〈*u*〉 must be of the same magnitude order inside the canopy [see also Wilson et al. (1998) for additional discussion]. In Fig. 4, the principal length scale computed from (4f), (18), and (20) is shown along with one–standard deviation bounds. Both length scale calculations (4f) and (18) are within the one–standard deviation bounds.

Another interesting finding is that the measured and Wilson’s (1988) model principal length scales (for *z* = 0.95*h*) are comparable to *h*/3 in the vicinity of *z*/*h* = 1.0. This vertical length scale is in excellent agreement with the length scale predicted by the mixing layer analogy of Raupach et al. (1996) for a “generic canopy.” Such length scale measurements suggest that vertical exchanges at the canopy–atmosphere interface do resemble a mixing layer. Nonetheless, our measurements suggest that the mixing layer analogy does not describe fully the horizontal length scale at *z*/*h* = 1.0. For example, anisotropy in eddy sizes, computed from *I*_{u}/*I*_{w} (representing *L*_{u}/*L*_{w}), is shown with bounds of one standard deviation, where *I*_{u} and *L*_{u} are integral timescales and length scales of the longitudinal velocity, respectively. Notice from Fig. 4 that the measured *L*_{u}/*L*_{w} is 10 and is half of its expected ASL value (e.g., Kader et al. 1989) at *z*/*h* = 1.0 but is still larger than would be expected by the arguments of Raupach et al. (1996). Inside the canopy, *L*_{u}/*L*_{w} decays to 3, and this estimate is more consistent with Raupach’s mixing layer analogy for a generic canopy.

In short, the vertical transport at *z*/*h* = 1.0 resembles a mixing layer but the horizontal transport does not. A possible explanation is attributed to the influence of inactive eddy motion (a boundary layer phenomenon) on the longitudinal velocity fluctuations (particularly the two-point statistics such as the correlation function or power spectrum). The influence of inactive eddy motion on the spectral characteristics of the longitudinal velocity (but not the vertical velocity) fluctuations at *z*/*h* = 1.0 was well demonstrated by Katul et al. (1997a, 1998) for the same forest stand. The role of inactive eddies is to increase *L*_{u}, as shown by the spectral analysis of Katul and Chu (1998).

## Conclusions

This study is the first to compare the performance of two higher-order canopy turbulence closure models for a forested canopy having a complex leaf area density profile. The comparisons demonstrated the following results.

The mean velocity profile is reproduced better by Wilson and Shaw’s (1977) model (about 45% improvement in rmse) when compared with Wilson’s (1988) model near the canopy–atmosphere interface.

The longitudinal velocity standard deviation profile is reproduced better by Wilson’s (1988) model (about 50% improvement in rsme) when compared with Wilson and Shaw’s (1977) model.

The modeled third moments by both closure schemes are inconsistent with the measurements. This failure is attributed to the flux gradient approximation used by both models to close the flux-transport term. Despite this failure, the Wilson (1988) model reproduced well all measured second-order statistics.

It is demonstrated that the longitudinal velocity standard deviation profile in the Wilson and Shaw (1977) model is sensitive to the parameterization of 〈

*w*′*u*′*u*′〉. The overestimation in computed 〈*u*′*u*′〉 by the Wilson and Shaw (1977) model is due to the poor parameterization of 〈*w*′*u*′*u*′〉. Similarly, the 〈*u*〉 overestimation by the Wilson (1988) model near*z*=*h*is due to the poor parameterization of 〈*w*′*u*′*w*′〉.The principal length scales of both Wilson’s (1988) and Wilson and Shaw’s (1977) models agree with integral length scale measurements inside the canopy. Both models and measurements result in a near-constant length scale in the middle third of the canopy where the leaf area density is highest.

While the measured vertical length scale

*L*_{w}and the length scale estimated from the Wilson (1988) model are consistent with the Raupach et al. (1996) mixing layer analogy, the measured horizontal length scale*L*_{u}is not. The anisotropy in length scales, measured by*L*_{u}/*L*_{w}, is larger by a factor of 3 than that suggested by Raupach et al. (1996) for a generic canopy. The overestimation in*L*_{u}relative to a mixing layer value is attributed, in part, to the role of inactive eddies near the forest–atmosphere interface.

## Acknowledgments

The authors would like to thank Cheng-I Hsieh for his help in the data collection and for his comments on the manuscript, John D. Wilson for the helpful discussion on the Wilson closure model, Roger Shaw for his comments and suggestions on the second-order closure model of Wilson and Shaw, and Betsy Nettleton for the many editorial suggestions. Part of this work was performed when the first author was on a semester sabbatical from the Nicholas School of the Environment, Duke University. This project was funded, in part, by the National Science Foundation (NSF-BIR 9512333) and by the Department of Energy through the FACE-FACTS project (DE-FG05-95 ER 62083). The FORTRAN-90 source codes for the two closure models can be obtained from the authors upon request.

## REFERENCES

_{2}and sensible heat densities and scalar profiles over and within a soybean canopy. Bound.-Layer Meteor., 61, 113–144.

_{2}in a pine forest under free air CO

_{2}enrichment. Oecologia, 104, 139–146.

_{2}sources, sinks, and fluxes in a uniform loblolly pine (Pinus taeda L.) stand. J. Geophys. Res., 102, 9309–9321.

_{2}concentrations to fluxes at local and regional scales. Aust. J. Bot., 40, 697–716.

### APPENDIX

#### Numerical Techniques and Closure Constants

##### The computational grid

The computational flow domain was set from zero to the highest measurement level (=1.44*h*). The grid node spacing is Δ*z* = 0.05 m, resulting in 415 grid nodes. This grid density was necessary because of rapid variability in leaf area density with distance.

##### Numerical scheme

The five ordinary differential equations (ODEs) for the Wilson and Shaw (1977) model [(4a)–(4g)] and the Wilson (1988) model [(14a)–(14e)] first were discretized by central differencing all derivatives. An implicit numerical scheme was constructed for each ODE with boundary conditions specified by (15). The tridiagonal system resulting from the implicit forms of each discretized equation was solved using the TRIDAG routine in Press et al. (1992, 42–43) to produce the turbulent statistic profile. Profiles for all variables were assumed initially, and a variant of the relaxation scheme described in 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, 107–111) to ensure finite second derivatives of *a*(*z*). Convergence was achieved when the maximum difference between two successive iterations in 〈*u*〉 did not exceed 0.0001%. We ran both closure models with Δ*z* = 0.02 and 0.1 m and checked that all solutions were independent of Δ*z.*

##### Closure constants

###### Wilson and Shaw’s (1977) model

For the Wilson and Shaw (1977) model, the closure constants are dependent on the choice of the boundary conditions and are determined by assuming that in the atmospheric surface layer (*z* > 2*h*), the flux-transport term is negligible and that 〈*u*′*w*′〉, 〈*u*′^{2}〉^{1/2}, and 〈*w*′^{2}〉^{1/2} become independent of *z* for near-neutral conditions. These simplifications result, after some algebraic manipulations, in the following relationships of *A*_{u}, *A*_{υ}, and *A*_{w} with *a*_{2}, *a*_{3}, and *C*:

where *A*_{q} = (*A*^{2}_{u} + *A*^{2}_{υ} + *A*^{2}_{w})^{1/2}. The closure constant *a*_{1} is determined by noting that the eddy diffusivity is *k*_{υ}(*z* − *d*)*u*∗ in the surface layer. Hence, *qλ*_{1} becomes identical to *k _{υ}*(

*z*−

*d*)

*u*∗ leading to

*a*

_{1}= 1/

*A*

_{q}. The above equations are the first analytic expressions that relate closure constants to ASL boundary conditions for the Wilson and Shaw (1977) model.

###### Wilson’s (1988) model

After the flux-transport terms are neglected; the wall function is noted to be unity in the ASL (i.e., *β* = *c*^{′}_{1}; the ASL formulations for 〈*u*′^{2}〉^{1/2}, 〈*υ*′^{2}〉^{1/2}, and 〈*w*′^{2}〉^{1/2} are replaced; and 〈*u*′*w*′〉 (=−*u*^{2}_{∗}) is assumed to be independent of height; the Wilson (1988) equations for 〈*u*′*w*′〉, 〈*u*′^{2}〉, 〈*υ*′^{2}〉, and 〈*w*′^{2}〉 reduce to

This system of equations is identical to the set of equations derived by Gibson and Launder (1978) for zero Richardson number [see their (26)–(29)]. The above set of equations is overspecified (i.e., it has more equations than unknowns) with a trivial solution *c*_{2} = 1, *c*^{′}_{1} = 0, and *c*_{1} = 0. This solution ensures a balance between mechanical production and dissipation but results in zero pressure redistribution. Gibson and Launder (1978) admit that, while their model cannot match the neutral atmospheric surface layer, it accounts for a wide range of flow types, including free shear and boundary layer flows. They determined *c*_{2}, *c*^{′}_{1}, and *c*_{1} optimally from laboratory measurements from both free-shear and boundary layer flows. It is important to note that their laboratory boundary layer measurements did not exhibit the same degree of anisotropy (as measured by *A*_{u}/*A*_{w}) as was seen in the near-neutral atmospheric boundary layer measurements reported in Wilson (1988) and this experiment.

Wilson (1988) used values for *c*_{2}, *c*^{′}_{1}, and *c*_{1} that were identical to those of Gibson and Launder (1978) even though his *A*_{u}, *A*_{υ}, and *A*_{w} were different from those in the Gibson and Launder study. Hence, for consistency with Wilson (1988) and Gibson and Launder (1978), the same values for the closure constants were used here. A summary of the resulting closure constants is presented in Table 1.

## Footnotes

*Corresponding author address:* Gabriel G. Katul, School of the Environment, Box 90328, Duke University, Durham, NC 27708-0328.