Large-eddy Simulations of Surface Roughness Parameter Sensitivity to Canopy-structure Characteristics

Surface roughness parameters, namely the roughness length and displacement height, are an integral input used to model surface fluxes. However, most models assume these parameters to be a fixed property of plant functional type and disregard the governing structural heterogeneity and dynamics. In this study, we use large-eddy simulations to explore , in silico, the effects of canopy-structure characteristics on surface roughness parameters. We performed a virtual experiment to test the sensitivity of resolved surface roughness to four axes of canopy structure: (1) leaf area index, (2) the vertical profile of leaf density, (3) canopy height, and (4) canopy gap fraction. We found roughness parameters to be highly variable, but uncovered positive relationships between displacement height and maximum canopy height, aerodynamic canopy height and maximum canopy height and leaf area index, and eddy-penetration depth and gap fraction. We also found negative relationships between aerodynamic canopy height and gap fraction, as well as between eddy-penetration depth and maximum canopy height and leaf area index. We generalized our model results into a virtual " biometric " parameterization that relates roughness length and displacement height to canopy height, leaf area index, and gap fraction. Using a decade of wind and canopy-structure observations in a site in Michigan, we tested the effectiveness of our model-driven biometric parameterization approach in predicting the friction velocity over heterogeneous and disturbed canopies. We compared the accuracy of these predictions with the friction-velocity predictions obtained from the common simple approximation related to canopy height, the values calculated with large-eddy simulations of the explicit canopy structure as measured by airborne and ground-based lidar, two other parameterization approaches that utilize varying canopy-structure inputs, and the annual and decadal means of the surface roughness parameters at the site from meteorological observations. We found that the classical representation of constant roughness parameters (in space and time) as a fraction of canopy height performed relatively well. Nonetheless, of the approaches we tested, most of the empirical approaches that incorporate seasonal and interannual variation of roughness length and displacement height as a function of the dynamics of canopy structure produced more precise and less biased estimates for friction velocity than models with temporally invariable parameters .


Introduction
Our ability to accurately predict mass and energy fluxes from the land surface to the atmosphere at any timescale depends on the accuracy of the surface drag parameterization (Finnigan, 2000;Mahrt, 2010).Over forested environments, vertical mixing of canopy air with the free atmosphere above, which is the process responsible for the exchange of energy, water vapor, and CO 2 between the land surface and the atmosphere, is a function of the turbulent eddies created through interactions between vegetative structure (e.g., trees, tree stems, leaves) and the wind (Thomas and Foken, 2007a).In many regional models, estimation of surface drag, and thus surface fluxes, is typically dependent upon parameterization of the friction velocity, u * , based on Monin-Obukhov similarity theory (MOST; Monin and Obukhov, 1954) using parameters that describe the effects of drag generated by the surface on the shape of the curve describing the ver-tical distribution of wind speed.These parameters are displacement height, d, and roughness length, z 0 .Though they represent different physical properties of the surface effects on the velocity profile, we will refer to them throughout the manuscript using the combined term "roughness parameters".In many land surface, vegetation, ecosystem, and hydrology models, such as the Community Earth System Model (CESM; Gent et al., 2011), Mapping Evapotranspiration with Internalized Calibration (METRIC; Allen et al., 2007), and Surface Energy Balance Algorithm for Land (SEBAL; Bastiaanssen et al., 1998), the surface sensible and latent heat fluxes are functions of the aerodynamic resistance for heat transfer, r ah .r ah is a function of the turbulence at the surface layer, defined through the friction velocity, u * .In models which cannot directly resolve u * , r ah is parameterized as a function of d and z 0 .In these models, d and z 0 may be derived from different canopy-structure characteristics.In the simplest approach, d and z 0 are linear functions of site-level canopy height (h) -typically d ≈ 0.66 h (Cowan, 1968) and z 0 ≈ 0.10 h (Tanner and Pelton, 1960).The accuracy of these estimates may be limited, however, by the dynamic nature (space and time) of canopy-structure characteristics.First, the canopy is a complex structure that is hard to describe using simple low-variable-number formulations.Second, estimates of the canopy-structure characteristics are limited by the typical absence of data about the vertical distribution of leaf area (Massman and Weil, 1999;Shaw and Pereira, 1982) and tree-top heights, as well as the difference between coarse model grid-cell resolution and the finer scale at which canopy-structure characteristics vary and affect roughness and momentum and flux transfer.
One common approach to incorporate canopy structure in the parameterization of roughness length into models in a more realistic way utilizes satellite imagery products to estimate vegetation structure and relate it to canopy-roughness relationships.For example, the SEBAL model (Moran, 1990) utilizes a function based on the normalized difference vegetation index (NDVI), while the METRIC model employs the Perrier function (Perrier, 1982).These canopy-roughness relationships have been shown to improve evapotranspiration estimates (Santos et al., 2012), but they are specific to sparse or short vegetative environments, such as agricultural systems, and are not typically recommended for forest environments (Bastiaanssen et al., 1998).
To incorporate the effects of canopy structure in denser and taller vegetative environments such as forests, empirical functions have been proposed using coarse canopy metrics such as canopy area index (the total, single-sided area of all canopy elements within a 1 m × 1 m ground area; Raupach, 1994), stand density (stems per area), or leaf area index (LAI, the total surface area of leaves found within a 1 m × 1 m vertical column of vegetation; Nakai et al., 2008a).However, the data required to use these functions are typically not available at most sites and, with the exception of LAI, are not yet obtainable through large-scale satellite remote sensing.In many climate models, surface-layer grid cells are prescribed with biome-specific qualities, i.e., sets of parameters describing constant vegetation structure and flux-driving characteristics for all model cells containing a specific biome or plant functional type (PFT).For example, the Ecosystem Demography model version 2 (ED2; Medvigy et al., 2009) provides 20 different vegetation functional types, 7 of which are representative of forested environments, to describe all land surfaces across the globe.Each such vegetation functional type is characterized by fixed, canopy-height-driven roughness parameters.Similarly, aerodynamic resistance to surface flux in the advanced hydrological model tRIBS+VEGGIE (Ivanov et al., 2008) is only driven by vegetation height, which is either prescribed or set as a default per PFT.
Roughness parameters have been shown to scale with structural characteristics, such as the influence of area-index (vegetation area per ground area) terms on d and z 0 , through numerical studies (Shaw and Pereira, 1982;Choudhury and Monteith, 1988) and wind-tunnel experiments (Raupach, 1994).Above-canopy meteorology data have shown estimates of roughness parameters to be highly variable both spatially and temporally (Maurer et al., 2013;Harman, 2012;Zhou et al., 2012).As evidence for canopy-roughness relationships has risen, various studies have attempted to generalize small-scale interactions between roughness parameters and canopy structure by deriving d and z 0 from above-canopy meteorological measurements (Braam et al., 2012;Maurer et al., 2013;Raupach et al., 1996;Nakai et al., 2008a), remote-sensing (Schaudt and Dickinson, 2000;Weligepolage et al., 2012), numerical experiments (Grimmond and Oke, 1999;Wouters et al., 2012), and large-eddy simulations (LESs; Aumond et al., 2013;Bohrer et al., 2009;Bou-Zeid et al., 2007, 2009).Although the understanding of these small-scale canopy-roughness interactions has grown, accounting for fine-scale canopy-structure effects on roughness parameters in larger-scale climate models requires further development.
In this study, we use the Regional Atmospheric Modeling System (RAMS)-based Forest Large-Eddy Simulation (RAFLES; Bohrer et al., 2008Bohrer et al., , 2009) ) to conduct a virtual experiment to estimate the sensitivity of surface roughness parameters to specific characteristics of fine-scale canopy structure.RAFLES incorporates a prescribed 3-D domain that includes the vegetation leaf density and stem diameters and dynamically calculates the change to wind velocity as a function of leaf and stem surface drag in each voxel (Chatziefstratiou et al., 2014).The level of detail at which vegetation is represented in RAFLES makes it particularly suitable for conducting this series of virtual experiments that simulate the drag parameters over a simplistic set of virtual canopy structures that vary by structural component, including stand density and patch fraction, canopy height, leaf area index, and vertical profile of leaf density.The approach of prescribing drag in LESs to resolve site-level roughness was previously tested and shown to provide higher accuracy than the tra-ditional roughness parameterization (Aumond et al., 2013).Finally, we use 10 years of direct observations of canopystructure and roughness parameters (Maurer et al., 2013) to estimate the sensitivity of modeled friction velocity to temporal variation in canopy structure and its effects on roughness length.We compare these results with other approaches that may be used to represent canopy structure when modeling roughness parameters.

Theory
Monin-Obukhov similarity theory (MOST) describes the relationships between the mean horizontal wind speed and the friction velocity in the inertial sublayer (Monin and Obukhov, 1954).In brief, MOST describes this relationship using a logarithmic function with parameters d and z 0 .Further details on the formulation of MOST used in this work are described in Maurer et al. (2013).The original MOST formulation was expanded to include the effects of thermal instability and the flow regime in the roughness sublayer (RSL), as follows: where u z is the mean horizontal wind speed at height z, above the ground.When the data are derived from meteorological observations, an overbar over a variable represents the 30 min mean of the 10 Hz time series of that variable.Given the mean eastward and northward wind velocities, u and v, u z is rotated toward the wind direction such that where κ is the von Kármán constant, ∼ 0.4, and z * is the upper limit of the RSL estimated as 2 h (Mölder et al., 1999;Raupach et al., 1996), with h the canopy height.I is an indicator function defined as (I = 1 for z ≤ z * ; or I = 0 for z > z * ).u * is the friction velocity defined as where each prime term (e.g., w ) is the perturbation of the specific variable from its mean (e.g., w − w).The atmospheric-stability correction function, ψ m (x), was described by Paulson (1970) for unstable atmospheric condi-tions (z/L < 0) as where x is either (z − d) /L or z 0 /L.Current understanding of aerodynamic properties near forest canopies within the roughness sublayer (RSL) has led to empirical corrections to the MOST model (Harman and Finnigan, 2007;De Ridder, 2010;Cellier and Brunet, 1992;Garratt, 1980;Mölder et al., 1999;Physick and Garratt, 1995;Raupach, 1992).These corrections allow us to utilize MOST with meteorological observation within the RSL, which typically includes the height range where eddycovariance measurements of forest flux dynamics are conducted across the globe.The RSL correction we used, ψ u (x 1 , x 2 ), was described by De Ridder (2010) as where ; and υ, µ, and γ are empirical constants provided by De Ridder (2010) as 0.5, 2.59, and 1.5, respectively.The inclusion of the RSL correction (ψ u = 0) occurs when the calculation is performed within the RSL (z ≤ z * , I = 1).Flux data are typically observed within the RSL at one point in space, requiring the implementation of the RSL correction.When boundary layer conditions are near neutral, (z − d) /L and z 0 /L approach zero, and thus ψ m (x) becomes negligible (Eq.4).Contrary to the classic estimate of z 0 (function of h), Thom (1971) suggested a relationship between z 0 and (h − d), as opposed to a relationship between z 0 and h alone, where the ratio of z 0 / (h − d) was defined as λ, a dimensionless, stand-specific parameter.This allows z 0 to be dependent on the spacing of the surface roughness elements and not only their height.For example, (h − d) will theoretically be smaller for more densely packed surfaces, providing a smoother surface and smaller roughness length.This relationship can be written as (6) Nakai et al. (2008b) substituted the aerodynamic height, h a , for the canopy height, h, into this relationship and rearranged the equation to read In simulation results, where the detailed 3-D wind field is known, we use Eq. ( 7) to calculate λ for each simulation using h a , which can be calculated from the vertical profile of horizontal wind speed and the empirically fitted d and z 0 .It is determined by identifying the height of the inflection point in the vertical wind-speed profile.This height marks the transition between the sub-canopy and above-canopy flow regimes (Thomas and Foken, 2007b).We investigated the eddy-penetration depth (δ e ), which is the length scale describing the vertical range from the top of the canopy that is influenced by turbulent mixing from above.It is defined as the distance between h a and the height where the momentum flux value is 10 % of its value at h a (Nepf et al., 2007).

Site description
The data used to test the effectivity of our LES-driven and other modeling approaches originate from a mixed, deciduous forest site at the University of Michigan Biological Station (UMBS) in northern, lower Michigan, USA (45 • 33 35 N, 84 • 42 48 W; elevation: 236 m above sea level).The forest is dominated (∼ 30 % of leaf area index) by early-successional bigtooth aspen (Populus grandidentata) and paper birch (Betula papyrifera), with a mean age of 85-90 years (Gough et al., 2013).The remaining leaf area is mostly represented by red oak (Quercus rubra), red maple (Acer rubrum), and white pine (Pinus strobus).Mean canopy height is roughly 20-25 m with an average stem density of ≈ 750 stems ha −1 (including only trees with diameter at breast height (DBH) > 8 cm).Eddy-covariance flux measurements have been ongoing at the site since 1999 and data are available through AmeriFlux (http://ameriflux.lbl.gov/);site code: US-UMB.Empirical allometric equations, fitted to measurements in this site (Garrity et al., 2012), are used to determine canopy height from a tree census and measurements of DBH.Full censuses were conducted in 2001 and 2010, and partial censuses of DBH for 993 trees are measured annually.Leaf area index is measured weekly using an optical sensor (LAI2000, LI-COR Biosciences, Lincoln, NE, USA).Additional details on the calculation of roughness length parameter from wind observations in the site and the determination of canopy structure are described in Maurer et al. (2013).Portable canopy lidar measurements (Hardiman et al., 2013) were used to determine the mean leaf area density profile that was used as the "natural" leaf area density case.Airborne lidar measurements were conducted by the National Center for Airborne Lidar Mapping (NCALM) in summer 2009.The lidar data and processing for our site are described in Garrity et al. (2012).This data set was used to determine the mean and variation of canopy top height and gap fraction, and to prescribe the explicit canopy structure in the "realistic" LES test case (see Sect. 2.4).

Large-eddy simulations
We used wind fields and heat fluxes from RAFLES simulations results to calculate surface roughness parameters of simplified virtual forests.RAFLES (Bohrer et al., 2009) uses a 3-D heterogeneous canopy domain where leaf and stem areas are prescribed within each voxel.The leaf area density and the instantaneous wind speed within the voxel determine the drag force that is applied to wind flow through that grid cell within each time step.Common to the approach used in most LESs, it assumes the leaf area is composed of flat surfaces oriented downstream and neglects higher-order effects of leaf and stem shapes and sub-grid-scale wake generation (shown to be a small effect; Shaw and Patton, 2003).It is combined with radiation attenuation (given the leaf densities in the grid cells above) to determine the sensible and latent heat fluxes emitted from each grid cell.The model uses the finite volume approach for discretization of the simulation domain.It resolves the effects of volume restriction due to the volume of the vegetation (stems, branches) by reducing the aperture areas available for flux exchange between each pair of neighboring grid cells and by reducing the volume that is available for flow within each grid cell according to the volume of the vegetation present (Chatziefstratiou et al., 2014).It resolves sub-grid-scale turbulence using the Deardorff (1978) scheme, and includes a parameterization for sub-grid-scale turbulence dissipation due to leaf drag (Shaw and Patton, 2003).
Simulations consisted of 3 h of simulation time at a time step of 0.02 s.RAFLES uses a nested time-stepping scheme with higher frequency calculations for turbulence and even higher frequency calculations for pressure perturbations.Eight pressure and four turbulence time steps were nested in one model time step.Output data snapshots of all grid cells in the simulation domain were recorded every 2 s.The initial 2.5 h of simulation time were used as a "spin-up" period to ensure satisfactory turbulent mixing and semi-stability of the vertical profiles of turbulence and potential temperature.The latter half hour of simulation time was used for analysis, consisting of 300 2 s snapshots.Synthetic virtual domains covered 1.25 km × 1.25 km × 1.4 km (length × width × height) at a horizontal grid spacing of 5 m × 5 m, which approximately corresponds to the mean size of individual tree crowns.Vertical grid spacing was 3 m in the lower subdomain, from the ground to 100 m above ground level.Above that region, vertical grid spacing was gradually increased by 12 % per each subsequent horizontal layer up to a maximal grid spacing of 30 m.The vertical grid spacing then remained constant above that height up to the model top at 1.4 km.The model has periodic boundary conditions at the lateral boundaries, no-slip boundary conditions at the bottom boundary, and a no-flux top boundary with Rayleigh friction to dampen vertical perturbations at the top six model layers (180 m).Initial conditions were horizontally homogeneous and followed a prescribed vertical profile for potential temperature, humidity, and wind speed.The prescribed initial vertical profile of the potential temperature described a well-mixed atmospheric boundary layer and was constant from 50 m to the height of the capping inversion, and increased with height above that level.Latent and sensible heat fluxes were prescribed based on observed mean noontime observations for August 2011 above the canopy at US-UMB.For each column of the horizontal simulation domain, the sum of the fluxes and Bowen ratio were distributed around the prescribed mean as an empirical function of LAI.Fluxes were further distributed vertically following a leaf-area-dependent empirical exponential profile.More details on the numerical setup of the model and the approach for flux forcing are provided in Bohrer et al. (2009).

Virtual experiment setup: sensitivity analysis to quantify the effects of specific canopy-structure characteristics on roughness parameters
Forest canopies are a complex array of 3-D structures.Many structural characteristics, such as tree height, LAI, vertical leaf area density (LAD) profile, and gap fraction, among others, affect the airflow inside and above the canopy and, consequently, the resulting roughness parameters and aerodynamic properties of the surface that describe such canopy structure.Using synthetic cases representing different aspects of canopy structure, we conducted a virtual experiment to test the sensitivity of roughness parameters to four axes of canopy structure: (1) mean site-level LAI, ranging from observed leaf-off conditions (LAI = 1.0 m 2 m −2 ) to typical, mid-growing season leaf-on conditions (LAI > 1.0 m 2 m −2 ); (2) LAD (m 2 m −3 ) profile, defined through the vertical bias of the vertical leaf density distribution (See Fig. S1 in the Supplement); (3) canopy height ranging from 9 to 27 m; and (4) canopy patch-level continuity (gap fraction) ranging from 0 to 50 % (see Fig. S2).Based on the available computing resources, we selected 20 combinations of the structural characteristics listed above.A list of all simulation cases and the canopy-structure characteristics is presented in Table 1.
In the gap fraction cases, canopy gaps were randomly created across the domain ranging from a single pixel (25 m 2 , tree-crown scale) to multi-pixel blocks (tens to hundreds of square meters).A gap was described by shorter vegetation (h = 9 m) and a non-gap (closed canopy) was described by taller vegetation (h = 27 m).It should be noted that we introduced gaps in our horizontally homogenous canopy using holes of varying sizes and shapes, which was done to minimize the complexity of the prescribed "heterogeneity" treatment (Fig. S2).The resulting gap-size distribution was arbitrary and may not have been well representative of an actual, heterogeneous canopy environment with tree-fall gaps.

Empirical determination of roughness parameters from simulations results
To calculate flux and wind statistics, we first calculated the mean value of each model variable at each vertical model level over the entire horizontal domain, and over all 300 time snapshots.We then rotated the horizontal wind coordinates of each vertical level toward the downstream direction, such that the resulting mean rotated downstream velocity is where xyt marks an average of the simulation results over all voxels in the x (eastward), y (northward), and t (temporal, 300 snapshots) dimensions.Although the wind forcing aloft is eastward, a rotation develops following the Ekman spiral and is further amplified by random x-y asymmetrices in the simulation domain.The rotation for the horizontal coordinate system of each horizontal layer is necessary to maintain a consistent downstream axis required for data analysis.After this rotation, we calculated the instantaneous perturbation of the velocity components from the xyt average for each voxel in space and time along each horizontal layer, such that where the prime indicates an instantaneous perturbation from the mean value in this example of the u r (downstream) velocity component.Similar formulation applies to the vertical (w) and cross-stream (v r ) velocity components.Momentum flux at the downstream direction was calculated as See Bohrer et al. (2009) for additional details on the calculation of wind statistics and momentum fluxes from RAFLES output.
We determined the effective aerodynamic canopy height, h a , by identifying the height of the inflection point in the vertical wind-speed profile as mentioned previously.To find this point, we compiled a domain-averaged wind-speed profile using Eq. ( 8).Then, we determined h a as the location where the second derivative of the horizontal wind profile crosses zero.We approximated this location within the vertical grid resolution using linear interpolation.We calculated the characteristic domain-averaged u * for each simulation case by calculating the horizontal-temporal average u * for each horizontal plane of grid cells within the 3-D virtual domain and further averaging these vertically over the range from 3.5 to 4.5 h (u * values are nearly invariable with height in that range).Obukhov length was calculated for each horizontal plane of grid cells within the 3-D virtual domain as a function of the characteristic u * , surface heat flux (prescribed), and the mean potential virtual temperature at each horizontal plane of grid cells.Next, the vertical profile of horizontal mean wind speed from all grid layers above 1.5 h a and below 4.5 h (95 m) above ground was fitted to Eq. (1) to determine d and z 0 using the characteristic friction velocity and the Obukhov length.The empirical fit was calculated using MATLAB's (version R2013b, The MathWorks Inc., Natick, MA, USA) nonlinear, least-squares fit function: fit.We constrained the solution for the surface roughness parameters to a physically meaningful range by constraining d to be between 0 and h a of the simulated forest and z 0 to be larger than 0.

Virtual experiment to explore canopy-roughness relationships
We found that d was significantly affected by maximum canopy height (h max ; three-way ANOVA, Table 2).We also found that h a and δ e were significantly affected by h max , LAI, and gap fraction (GF; Table 2).z 0 was not found to be significantly affected by any single aspect of canopy structure investigated within this study.As suggested by Thom (1971) and Nakai et al. (2008b), we checked the relationship between z 0 and (h ad) and found a significant relationship (r 2 = 0.72, P < 0.001).We found a positive relationship between d and h max (fit forced through [0.0], Fig. 1).
Surprisingly, canopy gaps showed little effect on d.A higher correlation existed between d and h max (r 2 = 0.78) than between d and mean canopy height (r 2 = 0.48) across the gap fraction sensitivity analysis.There was little change to d with increasing gap fraction, except for the scenario with 50 % gap fraction in the leaf-on simulations, which was significantly lower.Therefore, the relationship with h max (which was constant as the number of gaps increased) was selected instead of mean canopy height (which decreased as the number of gaps increased).Seasonality (leaf-on vs. leaf-off) also showed surprisingly small differences in d as height was varied, which had previously been observed at US-UMB (Maurer et al., 2013).We found positive h a − h max and h a − LAI relationships and a negative h a − gap fraction (GF) relationship (Fig. 2).We note that a positive h a − h relationship was previously observed at US-UMB using 12 years of meteorological data and tree-growth censuses (Maurer et al., 2013).By utilizing the suite of RAFLES simulations, we empirically calculated a single canopy − h a relationship as where a = 0.06 m, b = (−)0.69m, and c = (−)0.11m.We found a negative δ e − LAI relationship and positive δ e − h max and δ e − GF relationships (Fig. 3).As expected, we found δ e to be consistently higher during leaf-off periods compared to leaf-on periods at corresponding heights and gap fractions as wind was better able to penetrate the subcanopy.Increased LAI intensified the effect of gap fraction on δ e as the slope of the leaf-on fit line was larger than that of leaf-off periods.
Relationships were empirically determined using roughness parameters from each RAFLES simulation, except for those with "unnatural" vertical LAD profiles (i.e., the "upper", "middle", and "lower" LAD cases) as no patterns were observed between any roughness parameters and vertical LAD profile.Maximum canopy height was used instead of mean canopy height because maximum canopy height was more tightly correlated with each roughness parameter than mean canopy height.The resulting roughness parameters for each simulation are listed in Table 1.
We calculated a "biometric" h a_b for the US-UMB site using the relationship we found in the virtual experiment between h a_b and LAI, gap fraction, and h max (Eq.12).To simulate the conditions in our site at US-UMB, we assumed a gap fraction of 5 %, which was found by calculating the percent area within the NCALM lidar scan domain with vegetation height less than 2 m.We used the peak growing season site-level mean LAI of 4.2 as measured from 2000-2011 (Maurer et al., 2013).A biometric d was then calculated using Eq. ( 11).Finally, a biometric z 0 was calculated as where λ = 0.34 was determined from Eq. ( 7) given the set of h a , d, and z 0 values from our simulations through the virtual sensitivity experiment.

Testing empirical approaches that link roughness parameters to biometric measurements
The biometric approach, derived from our simulation results, provides relationships between easily measurable characteristics of the canopy (i. and d and z 0 .In order to evaluate the potential improvement to estimates of u * using this approach, we compared the accuracy and precision of modeled u * values using the biometric approach with those of five alternative approaches.We evaluate the resulting friction velocities predicted by each of these six (biometric and five alternatives) structure-driven parameterization approaches using 30 min observed values of u * , canopy height, and LAI over multiple years at US-UMB (2000-2011, at 34 m a.g.l).The five alternative approaches employed are as follows: 1. "Classical" -fixed d = 0.66 h and z 0 = 0.10 h, where we use h = 22 m based on a long-term average over the flux footprint during observation period.
2. "Explicit-LES" -fixed d = 0.67 h and z 0 = 0.094 h as determined from the simulation results of the "realistic" LES case.

"
Yearly observed" -a purely empirical approach, using the values of d and z 0 calculated from meteorological observations during each growing season at US-UMB from 2000-2011 (Maurer et al., 2013).In this approach, the values of d and z 0 vary each year according to observations.d and z 0 were calculated by fitting Eq. 1 to a seasonal set of half-hourly mean observations of wind speed and friction velocity at twice the canopy height (46 m a.g.l.) and only during neutral to slightly unstable atmospheric conditions during daytime.We also tested applicability of shorter-term observations of d and z 0 to long-term predictions of friction velocity.This test was motivated by the fact that there are only a few sites around the world with more than a decade of data, while short observation campaigns are more common.
We used the observed d and z 0 from each year to simulate the entire decadal time series of friction velocity.This resulted in 12 different "yearly" models.Anecdotally, the most accurate model was associate with observed d and z 0 from 2008, and the least accurate model with the yearly values from 2005.
Numerous past studies have attempted to derive relationships between roughness parameters and other canopy-structure statistics.We chose two in this study: 4. Raupach (1994) calculated d and z 0 as functions of canopy area index ( ), drag coefficient (c d ), and canopy height (h): and where c d = 7.5, η h = 0.193, and = 2nbh/A, where n is the number of stems in a sample plot, b is the mean diameter at breast height, h is the mean tree height, and A is the total ground area within the canopy sampling area.Full plot censuses provided the data to calculate .These were conducted in 2001 and 2010, and values were linearly interpolated for the years between the censuses and extrapolated to 2011. 5. Nakai et al. (2008a) calculated d and z 0 as functions of stand density (ρ s ), LAI, and h: and where α and β are 7.24 × 10 −4 ha stems −1 and 0.273, respectively, and we used the US-UMB mean stand density of 750 stems ha −1 .
The values of d and z 0 as determined by each of the parameterization approaches are listed in Table 3.The range for yearly observed mean d values was 18.3-26.0m and for z 0 0.99-1.99m.The classical approximation based on h resulted in a significantly lower d = 14.0 m (outside the range of the interannual variability over 12 years), and a slightly above-range z 0 = 2.10 m.The explicit-LES approach resulted in a very similar d to the classical approach.The biometric approach predicted high but within-range d values (24.0-25.0m) but extreme z 0 values (3.64-3.82m).There was nearly no overlap between the values of z 0 from each of the approaches, indicating poor agreement between approaches for this parameter.
Table 3. Thirty-minute block-averaged friction velocity (u * ) model evaluation against measured u * for displacement height (d), and roughness length (z 0 ) calculated from various methods -"classical", "yearly observed", "biometric", "Raupach 94", and "Nakai 08" -at US-UMB spanning the 2000-2011 growing seasons.We show the slope and intercept of the linear fit, which are measures of the accuracy of the models; the coefficient of determination (r 2 ), which is a measure of precision; and the root-mean-square error (RMSE) between modeled and observed u * , which is indicative of both precision and accuracy.

Improvements to estimates of friction velocity using canopy-structure-roughness relationships
Modeled u * from all six approaches was regressed against observed u * .The slope and intercept of the fit line (estimates of accuracy), coefficient of determination (r 2 ), and root-mean-square error (estimates of precision) are reported in Table 3. Surprisingly, all parameterization approaches produced similar results, with coefficients of determination between 0.56 and 0.61, near zero, but significantly negative intercepts between (−)0.052 and (−)0.072(significant margin ±0.004).The most significant difference between the approaches was in their bias.All approaches (except the yearly observed one of 2008, which was the only one that was not significantly biased) produced a significant positive bias, but the bias varied from near zero to 43 % (slope of observed vs. modeled fit line between 1.01 and 1.431, significant margin ±0.01).The results of all parameterization approaches are listed in Table 3.We found that the precision of the results obtained by using each of the 12 yearly observed models over the entire 12-year period to be higher than the combined results of using the observation for each specific year during that year only.The bias of the prediction obtained with the observed d and z 0 , applied to the entire 12-year period, varied from no significant bias (using the 2008 parameters) to 1.38 (with the 2005 parameters).The combined (each year with its own parameters) produced an intermediate bias for the friction velocity estimates.
The yearly observed method is dependent on long-term observations of wind, temperature, heat flux, and friction velocity, which are rarely available in forest sites.The other methods we tested do not require directly observed roughness parameters.Of these methods, the "Raupach 94" approach had the highest precision and lowest bias (slope = 1.24, r 2 = 0.604), the explicit-LES approach ranked second and our biometric approach ranked third, although it performed similarly to the very simple classical approach.The "Nakai 08" approach proved to be the least compatible with our site.

Response of roughness parameters to canopy-structure change
To date, despite a strong need by the modeling community, there is no single consensus approach that relates roughness length and displacement height to observable properties of canopy structure, such as LAI, height, leaf density, and gap fraction.Furthermore, observations in our field site (Maurer et al., 2013) and by others (Nakai et al., 2008a) have shown that the roughness parameters in forests are not easily constrained by leaf area or canopy height.Our underlying assumption in setting up this model-based experiment was that the lack of a clear empirical relationship between roughness parameters and canopy structure was due to the complexity of canopy structure.We assumed that different characteristics of the canopy drive different effects on roughness length and displacement height.In real forests, many of the structural characteristics vary in time in different ways, resulting in interacting and sometimes conflicting effects on roughness length and displacement height.We set up a numerical experiment that was designed to separate the effects of different observable characteristics of canopy structure.We also hypothesized that, to some degree, the difficulty in identifying a clear effect of canopy structure on each of the roughness parameters is because roughness length and displacement height values may trade off, such that similar solutions can be fitted either with low d or high z 0 , or vice versa (Nakai et al., 2008a, b;Maurer et al., 2013).By testing the independent effects of different characteristics of canopy structure through a set of controlled virtual experiments, we indeed found that different roughness parameters where sensitive to different structural characteristics.The aerodynamic canopy height (h a ) and eddy-penetration depth (δ e ) were both sensitive to leaf area, canopy height, and gap fraction (Figs. 2, 3).In contrast, d was only significantly sensitive to canopy height, while z 0 did not show any significant relationships with any single canopy-structure characteristic.
We found positive d −h max and h a −h max relationships independent of LAI.A strong correlation had previously been reported between h a and h (Nakai et al., 2008b;Bohrer et al., 2009;Maurer et al., 2013;Thomas and Foken, 2007b).As canopy height was the only canopy characteristic that varied among the "canopy height variation" simulations (Table 1c), it is reasonable to assume that δ e would be relatively constant, regardless of canopy height.However, as canopy height increased within our virtual domain, the constant mean site-level LAI was stretched further in the vertical direction.Therefore, the mean leaf density in the upper canopy was smaller for taller canopies, resulting in an increased δ e with canopy height (Fig. 3b).In spite of increased δ e , we also observed a positive d −h max relationship, indicating that the increased δ e only partially compensated for the increase in canopy height, allowing for d to increase linearly with canopy height, but with a slope smaller than 1.
We found a linear relationship between h a and gap fraction.Eddy-penetration depth scaled with gap fraction as well.It was consistently larger during leaf-off periods compared to leaf-on periods, and the presence of higher LAI during the leaf-on periods resulted in a steeper linear slope of the relationship between δ e and gap fraction (Fig. 3c).Intuitively, increased gap fraction should lead to increased δ e , as more canopy openings allow eddies to penetrate deeper into the canopy.These findings are not surprising, as Shaw et al. (1988) found deeper δ e at lower LAI.For example, we found that increased LAI and increased gap fraction corresponded to increased horizontal wind speed, momentum flux, and turbulence inside the canopy, below 1 h (Figs. 4,5,6).This was likely due to the extension of turbulent eddy penetration deep into canopy gaps, indicated by elevated standard deviation of the vertical velocity, σ w (a component of the turbulence kinetic energy), in canopy gaps (Fig. 6a).Such locations of increased turbulent eddy penetration are less likely to occur in horizontally homogenous canopies (Fig. 6b).However, the lack of any relationships between roughness length and gap fraction at all levels below 50 % gap (Table 1) was surprising, as Bohrer et al. (2009) found increases to d, z 0 , and h a in patchier canopies (more gaps) during leaf-on conditions.The major difference between these two studies was that the scale of the gaps prescribed here -corresponding to 1-2 crown sizes -was typically smaller than those in the Bohrer et al. (2009) experiments.We found no consistent correlations between roughness parameters and the mode of the vertical LAD profile, as the variability in roughness parameters over the range of LAD scenarios was extremely high (Table 1).Although the shape of the vertical profile of wind speed is apparently different between the "lower" and the "upper" LAD profiles (Fig. 7) there was no consistent canopy-wind or canopy-turbulence relationships that could be predicted by the bias of the vertical LAD curve (Fig. 7).LAD profiles may change in complex ways across the landscape and over many timescales (seasons, years, decades) due to disturbance or senescence.As our virtual experiment has shown, the effects of the vertical LAD profile are inconsistent with a simple representation of the vertical distribution of LAD using its vertical bias as a single descriptive characteristic.Our results indicate that site-level mean LAI and canopy height are easier to obtain and, in general, provide more reliable characteristics of canopy structure than the vertical profile of LAD.
Our simulations did not detect a continuous increase to d or z 0 with LAI, which was inconsistent with several previous wind tunnel or model studies (Choudhury and Monteith, 1988;Grimmond and Oke, 1999;Raupach, 1994;Shaw and Pereira, 1982).We also did not find significant relation- ships with any single property of canopy structure, except between displacement height and canopy height.To a limited degree, this was the result of tradeoffs between the two, as indicated by the fact that h a , which combines d and z 0 through the slope of their tradeoff curve, λ, was better constrained than d or z 0 alone.However, this tradeoff cannot fully explain the lack of relationship, as we did not find a significant and consistent relationship between z 0 and different canopy-structure characteristics even when we assumed a fixed displacement height and fitted only for z 0 (results not shown).Combined, our results indicate that both of our underlying hypotheses were at least partially false, and neither the structural complexity of the canopy nor the tradeoffs between z 0 and d can fully explain the lack of clear relationship between canopy structure and d and z 0 .
The lack of canopy-structure effects on z 0 within the virtual sensitivity experiment and, in particular, the lack of consistent seasonal differences between leaf-on and leaf-off periods may suggest that leaf area is not the primary driver of z 0 .To further understand the drivers of z 0 , we calculated the sensitivity of z 0 to changes in wind speed at a measurement height z above the canopy, δ z0u .This can be done by solving Eq. (1) for z 0 assuming neutral conditions, and calculating the sensitivity as the partial derivative of z 0 with respect to u z : We determine that, at low to intermediate mean wind speeds (below 3 m s −1 ), z 0 is extremely sensitive to variation in u, with the derivative being between 5 and 30 (Fig. 8).This indicates that, for an observed variation of 0.1 m s −1 measured at twice the canopy height, the resulting z 0 will change by 0.5-3 m, which is a full range of the expected z 0 values for a 20 m tall canopy.At our site in Michigan, 3 m s −1 was approximately the median wind speed and was therefore selected to drive the simulations.In reality, variations in halfhourly mean wind speed on the order of 0.1 m s −1 can be a result of local variations in the flow field due to topography or measurement errors due to instrument placement and calibration.In both reality and LES, such variations in wind speed at a given measurement point could also be the result of effects of local modification to the flow field due to specific heterogeneous canopy-surface structures (which were determined to extend up to 5 h; Raupach and Thom, 1981;Bohrer et al., 2009), and could also be driven by random large eddies that may affect the 30 min average at a specific half hour.We hypothesize that this high sensitivity of z 0 may be inhibiting the attempts to empirically estimate its relationships with the canopy-structure characteristic.

Integrating canopy-structure characteristics into models
Typically, surface roughness parameterization is used in models to directly or indirectly predict the friction velocity, which is further used in the surface flux calculations.
To test the performance of different parameterization approaches, we used data from 12 years of wind, friction ve- locity, Obukhov length, and canopy-structure observations in a forest site in Michigan.We compared six approaches that differ in whether (or not) they incorporate temporal variation into canopy structure, and in the source of data they require to determine z 0 and d.Surprisingly, but optimistically for the purpose of accurate modeling, all the surface roughness parameterization approaches we tested resulted in relatively high precision (r 2 = 0.58-0.61) in predicting the half-hourly friction velocity over 12 years.This is surprising because each of the approaches used a different set of values for z 0 and d, which in some cases were very far from each other.For example, the biometric and the classical approach performed rather similarly, but the biometric approach z 0 values were about 80 % larger than those of the classical approach.To understand this discrepancy, we calculated the sensitivity of the friction velocity to variation in z 0 , δ u * z0 .
For a case similar to the one we simulated, with a canopy at 22 m and mean wind speed of 3 m s −1 , we found that the friction velocity is not sensitive to changes in roughness length when roughness length is higher than 0.6 m (Fig. 8).As a general approximation (following the classical approach), for a forest canopy higher than 10 m, roughness length is expected to be larger than 0.1 h = 1 m.Therefore, while the value of the roughness length parameter is highly sensitive to changes in the half-hourly mean wind speed (Eq.17, Fig. 8, Table 1), the resulting friction velocity may not be greatly affected by this variation in the parameter's value.).We plotted the response curve over the same parametric range and the expected range of values for z 0 , when mean wind speed is around 3 m s −1 .u * is relatively insensitive (δ u * z0 < 0.15) for any z 0 above 0.5 m.
The best performing approach for parameterization of roughness length and displacement height was obtained using the annually observed values of these parameters.The yearly observed model demonstrated ∼ 7 % less error than the fixed-in-time classical canopy-roughness relationships.The combined yearly observed approach used the z 0 and d values for each year to predict friction velocity values in the same years.This method performed better than when applying the data observed during a single year to the entire time period.However, the roughness parameters observed during 2011 provided a more accurate and precise model for the entire 12-year time series than the combined approach.The z 0 and d values observed during 2005 provided the worst model, but they still performed better than the classical approach.It is rather intuitive that, when observations of z 0 and d exist, they will provide the best approach for modeling of friction velocity (Table 3).Our results indicate that the interannual variability in canopy structure that affects roughness length has only a very small effect on the resulting friction velocity.Annual growing-season averages of z 0 and d from any single year can provide a suitable approximation to the decadal time series of roughness length parameter values.However, the low spatial coverage by flux networks over the globe limits the use of this method across large spatial domains.
LESs with an explicit, prescribed canopy structure based on lidar observations of the canopy at a site can a surrogate virtual observations from which to evaluate the roughness parameters.However, these types of simulations are limited in their temporal domain (just a few hours as a representative of an entire decade).They are also dependent on high-resolution canopy lidar observations, which, to date, are not common.Parameterization approaches which rely on biometric observations, rather than on wind observations, may be the most reliable and broadly available method to estimate long-term roughness parameters.Our ability to estimate canopy-structure characteristics such as LAI, canopy height, and gap fraction over a broad range of spatial and temporal scales is continuously improving through the use of on-site biometric measurements as well as airborne and satellite remote sensing observations (Chen et al., 2002;Jonckheere et al., 2004;Zheng and Moskal, 2009).
As an indication for the potential of biometric approaches, the approach suggested by Raupach (1994) performed even better than the combined yearly observed approach (Table 3).However, this approach relies on stem census observations.While such records are more common than flux sites, there is still no broad global coverage for this type of observation.We tested two biometric approaches that only required more commonly observable canopy characteristics.LAI, canopy height, and gap fraction or stand density are required by both the Nakai et al. (2008a) approach and the approach derived by the virtual experiments in this study (the biometric approach) in order to determine z 0 and d.Of the two, our biometric approach performed better, and also provided slightly better estimates than the classical approach.Variable success by the three biometric methods may not be surprising -a study by Grimmond and Oke (1999) determined that careful consideration must be given to higher-order structural features of the surface than the ones represented in this study and include in the biometric approaches.Examples of such higher-order structural characteristics include the complexity of organization and the density of roughness elements.Similar reasoning could provide insight regarding the poor performance of the method of Nakai et al. (2008a) at US-UMB, which is less dense, taller, and has higher LAI than those sites used to parameterize the Nakai 08 method.
The biometric method presented in this study is essentially a variant of the classical method, with the major difference being the use of a variable maximum canopy height as opposed to mean canopy height, and adding small perturbations to roughness length based on LAI and gap fraction .The limited success of this method can be attributed to some degree to the limited effect of interannual variability in canopy structure.However, a decade of observations in a site represents only a very narrow range of potential canopy structures.Our simulation results suggest that this method could potentially improve the prediction of friction velocity when applied to situations where canopy-structure variability is larger, such as after significant disturbance events.

Conclusions
In this study we used an LES, long-term meteorological observations, and remote sensing of the canopy to explore the effects of canopy structure on surface roughness parameters in a forest site.We performed a virtual experiment to test the sensitivity of roughness parameters with respect to four axes of variation in canopy structure: (1) leaf area index, (2) the mode of the vertical profile of LAD, (3) canopy height, and (4) gap fraction.We found consistent relationships between the aerodynamic canopy height and LAI, maximum height, and gap fraction, as well as between d and maximal canopy height.We found that the predicted values of friction velocity are not sensitive to roughness length.As a result, most of the roughness-based approaches we tested for simulating friction velocity performed similarly well.This is despite having very different approaches for determining the values of z 0 and d, and having large differences in the range of z 0 and d values.This is good news for modelers, because it limits the error from using the current approaches that do not vary in time and do not incorporate canopy structure.
Nonetheless, most of the approaches we tested which used annually variable z 0 and d and that incorporated canopy structure provided better approximation for friction velocity than the classical, time-invariable method.Many easily obtainable metrics of canopy-structure characteristics are available through a suite of measurements, such as on-site meteorological and biometric observations or satellite-derived site characteristics.Additionally, many ecosystem models and ecosystem modules within Earth system models resolve the growth of the forest and accurately predict canopy height and LAI.Some models, such as the Ecosystem Demography model (Medvigy et al., 2009), even resolve the distribution of stem sizes.Such demographic models could readily incorporate the approach by Raupach (1994) for a significant improvement in surface roughness parameterization.For other models that resolve, or are forced by observed leaf area and vegetation height, our LES-derived biometric approach could offer an easy way to dynamically affect the roughnesslength parameterization.This could provide an improvement of surface flux modeling, especially when canopy-structure variations are large.Due to limited spatial coverage by direct meteorological measurements, remotely sensed structure statistics, and stand inventories, we suggest utilizing site-and time-specific biometric measurements of canopy structure to estimate site-level d and z 0 .The effectivity of these model improvements will, of course, be dependent upon the quality, quantity, and resolution of the data sets available at the forest of interest.
The Supplement related to this article is available online at doi:10.5194/bg-12-2533-2015-supplement.

Figure 5 .
Figure 5. Vertical profiles of (a) horizontal wind normalized by friction velocity and (b) momentum flux normalized by the square of friction velocity in a 27 m tall canopy with gap fractions of 0 % (blue), 10 % (cyan), 25 % (green), 35 % (orange), and 50 % (red), as well as in a continuous 21 m tall canopy (dashed back).Canopy height for the tall and short canopies is shown as dashed horizontal gray and green lines, respectively.

Figure 6 .
Figure 6.Vertical cross section through the simulation results of (a) a 27 m tall canopy with 25 % gap fraction and (b) homogeneous 21 m tall canopy.Thirty-minute mean wind speed and direction are illustrated using black arrows, the standard deviation of vertical velocity (an indication of turbulence intensity) is plotted using a color map.Canopy top in each simulation is illustrated by a solid green line.

Figure 7 .
Figure 7. Vertical profiles of (a) horizontal wind normalized by friction velocity and (b) momentum flux normalized by the square of friction velocity for "lower" (blue), "middle" (cyan), "upper" (green), and 'natural" (orange) LAD profiles.Canopy height shown as dashed horizontal green line.
Figure 8. (a)Sensitivity analysis of z 0 as a function of variation in the mean wind speed (δ zou ).We illustrate it here over a particular range of parameters, choosing a canopy height of h = 22 m (roughly the height we used in the simulation and observation site), displacement height of d = 0.67 h, observation height of 2 h (the recommended observation height for a flux tower), and u * of 0.35 m s −1 .The results are similar for other canopy heights and u * values.(b) Sensitivity of u * to variation in z 0 (δ u * z0 ).We plotted the response curve over the same parametric range and the expected range of values for z 0 , when mean wind speed is around 3 m s −1 .u * is relatively insensitive (δ u * z0 < 0.15) for any z 0 above 0.5 m.

Table 1 .
Description of simulation cases used for sensitivity analysis of roughness parameters derived from an LES over variable canopy layouts, and the resulting roughness parameters for each simulation case.Canopy structure was varied along four axes -(a) LAI, (b) vertical LAD profile, (c) canopy height, and (d) gap fraction -and included an additional (e) realistic simulation case.

Table 2 .
Results of a three-way ANOVA to test any significance that maximum canopy height (h max ), leaf area index (LAI), and gap fraction (GF) have on displacement height (d), roughness length (z 0 ), aerodynamic canopy height (h a ), or eddy-penetration depth (δ e ).P values listed in bold font indicate a significant effect.