1. Introduction
A forest is defined as an area of more than 0.5 ha, of which over 
$10\,\%$ is covered with trees with a minimum height of 5 m. Still according to the Global Forest Resource Assessment (FAO 2012, 2016), forests cover 
$30.6\,\%$ of Earth's total land area and their biomass and soils store more carbon than the atmosphere (Pan et al. Reference Pan, Birdsey, Phillips and Jackson2013). It is not surprising, therefore, that the flow within forests has been a subject of research, and the review papers by Raupach & Thom (Reference Raupach and Thom1981), Finnigan (Reference Finnigan2000), Belcher, Harman & Finnigan (Reference Belcher, Harman and Finnigan2012) and Brunet (Reference Brunet2020) offer an overview of how our understanding and modelling of forest flows have developed.
Initial studies (cf. Raupach, Thom & Edwards Reference Raupach, Thom and Edwards1980) considered a forest to be a perturbation to the upstream boundary layer. This approach (based on the characteristic roughness, 
$z_0$, and displacement height, 
$d$) can only be successfully applied to the mean flow over large homogeneously forested regions and is most useful under conditions of insufficient spatial resolution (e.g. Silva Lopes, Palma & Piomelli Reference Silva Lopes, Palma and Piomelli2015). Nowadays, the flow over a forest is viewed as an analogue of a mixing layer (cf. Finnigan Reference Finnigan2000), composed of large turbulent structures that gradient-diffusion-based approaches cannot mimic. Within the forest, the flow between the complex network of trunks and leaves increases shear and, consequently, turbulence – wake production – while, at the same time, large eddies are broken into smaller ones, increasing dissipation – spectral shortcut (e.g. Novak et al. Reference Novak, Warland, Orchansky, Ketler and Green2000; Cava & Katul Reference Cava and Katul2008).
In the case of a higher resolution, trees are seen as a drag force (
$f_i$), and modelled via additional contributions to the transport equations for momentum (
$U_i$, (1.1)), turbulent kinetic energy (TKE; 
$k$, (1.2)) and its dissipation rate (
$\varepsilon$, (1.3)). Traditionally, the canopy models in the Reynolds-averaged Navier–Stokes (RANS) formulation have the general form
\begin{gather} \mathcal{F}^{k-\varepsilon}_{U_i} = - C_d a(z) | \boldsymbol{U} | U_i, \end{gather}
\begin{gather}\mathcal{F}_k^{k-\varepsilon} = C_d a(z) \left( \beta_p | \boldsymbol{U} |^3 - \beta_d | \boldsymbol{U} | k \right)\!, \end{gather}
\begin{gather}\mathcal{F}_{\varepsilon}^{k-\varepsilon} = C_d \, a(z) \left( C_{\varepsilon 4} \beta_p | \boldsymbol{U} |^3 \frac{\varepsilon}{k} -C_{\varepsilon 5} \beta_d | \boldsymbol{U} | \varepsilon \right)\!. \end{gather}
Here 
$C_d$ is a mean drag coefficient, 
$a(z)$ is the leaf area density (
$LAD$), a function of the distance to the ground, and 
$|\boldsymbol {U}|$ is the mean scalar velocity, 
$( U_i U_i )^{1/2}$. Here, we are assuming the mathematical formulation of canopy flows in the context of the two-equation (
$k$–
$\varepsilon$) RANS model, first presented by Svensson & Häggkvist (Reference Svensson and Häggkvist1990) and Green (Reference Green1992). Since this model is based on the eddy-viscosity concept, its accuracy, when compared with Reynolds-stress models, is questionable in strongly advective conditions. However, it is successfully used in the study of atmospheric flows and by the wind industry (e.g. Costa et al. Reference Costa, Castro, Palma and Stuart2006; Dalpé & Masson Reference Dalpé and Masson2008; Abiven, Palma & Brady Reference Abiven, Palma and Brady2011).
Equations (1.2) and (1.3) were obtained using dimensional relations and the coefficients determined as usual in turbulence modelling; i.e. on the basis of canonical flows, with comparison of the computational and experimental results of mean wind speed, shear stress or TKE (see, for instance, Wilson Reference Wilson1988; Lien & Yee Reference Lien and Yee2004; Lien, Yee & Wilson Reference Lien, Yee and Wilson2005; Sanz & Katul Reference Sanz and Katul2007; Sogachev et al. Reference Sogachev, Kelly and Leclerc2012) or of large-eddy simulations (LES) (Silva Lopes et al. Reference Silva Lopes, Palma and Viana Lopes2013). Experimental data within real forests are scarce, and many of the constants have been fine-tuned with measurements taken outside forests, where the impact of such a different model (or set of constants) is lower. These equations try to accommodate the two opposing and competing phenomena of wake production and spectral shortcut. This effort, eloquently summarized in table 1, has led us to the current position, in which, in the context of two-equation turbulence modelling, neither standard parametrization nor standard values for the coefficients are available when modelling the flow over forests. Sanz (Reference Sanz2003), Sogachev & Panferov (Reference Sogachev and Panferov2006), Sanz & Katul (Reference Sanz and Katul2007), Sogachev (Reference Sogachev2009) and Sogachev et al. (Reference Sogachev, Kelly and Leclerc2012) provide some of the most comprehensive derivations of the coefficients, based on such phenomenological arguments as assuming the Kolmogorov relation, a dense and homogeneous canopy and constant mixing length within the canopy. Nevertheless, these physical constraints are insufficient and the coefficient values are found to be bound to such additional ad hoc restrictions as 
$C_{\varepsilon 4}=C_{\varepsilon 5}$ (Sanz Reference Sanz2003).
Table 1. Models and model coefficients in RANS formulation. Optimal values are underlined.
$^{a}$Based on Wilson et al. (Reference Wilson, Finnigan and Raupach1998).
$^{b}$Based on Svensson & Häggkvist (Reference Svensson and Häggkvist1990).
$^{c}$Based on Yamada (Reference Yamada1982) and Svensson & Häggkvist (Reference Svensson and Häggkvist1990).
$^{d}$Based on Green (Reference Green1992) and Liu et al. (Reference Liu, Chen, Black and Novak1996).
$^{e}$Based on Sogachev & Panferov (Reference Sogachev and Panferov2006).
$^{f}$Based on Sanz (Reference Sanz2003).
$^{g}$The mathematical formulation by Wilson (Reference Wilson1988) differs from the other models.
The forest flow parametrizations in the LES and RANS frameworks are highly similar. Both require a canopy-related term in the momentum equation and in the TKE (
$e$) of the subgrid scales (cf. Sagaut Reference Sagaut2006), if such an equation (1.4) is required (e.g. Shaw & Schumann Reference Shaw and Schumann1992). According to Kanda & Hino (Reference Kanda and Hino1994),
\begin{equation} \mathcal{F}_e^{LES} = C_d a(z) \left(\alpha |\boldsymbol{U}|^3 -2|\boldsymbol{U}| e \right)\!, \end{equation}
where 
$\alpha$, the ratio for the transformation of the energy loss in grid-scale flow through leaf drag to wake production in the subgrid-scale flow, can vary between 0 and 1. This equation is analogous to (1.2) and the coefficients 
$\alpha$ and 2 are denoted by 
$\beta _p^{LES}$ and 
$\beta _d^{LES}$. It should be mentioned that this analogy between (1.2) and (1.4) is merely formal, since 
$k$ encompasses the whole range of scales whereas 
$e$ is the TKE of subgrid scales. Most LES-based studies have their origin in the work of Shaw & Schumann (Reference Shaw and Schumann1992) and have used a simpler form of (1.4), in which 
$\beta _p^{LES}$ is equal to 0 and 
$\beta _d^{LES}$ is equal to 8/3 (e.g. Shaw & Patton Reference Shaw and Patton2003; Bohrer et al. Reference Bohrer, Katul, Walko and Avissar2009; Chatziefstratiou, Velissariou & Bohrer Reference Chatziefstratiou, Velissariou and Bohrer2014; Nebenführ & Davidson Reference Nebenführ and Davidson2015) or 2 (e.g. Shaw & Schumann Reference Shaw and Schumann1992; Kanda & Hino Reference Kanda and Hino1994; Dwyer, Patton & Shaw Reference Dwyer, Patton and Shaw1997; Shen & Leclerc Reference Shen and Leclerc1997; Patton et al. Reference Patton, Shaw, Judd and Raupach1998; Yang et al. Reference Yang, Morse, Shaw and Paw U2006; Dupont & Brunet Reference Dupont and Brunet2008, Reference Dupont and Brunet2009; Schrottle & Dornbrack Reference Schrottle and Dornbrack2013; Kanani et al. Reference Kanani, Träumner, Ruck and Raasch2014; Schlegel et al. Reference Schlegel, Stiller, Bienert, Maas, Queck and Bernhofer2015; Boudreault et al. Reference Boudreault, Dupont, Bechmann and Dellwik2017). Unlike RANS-based studies (see table 1), very few LES-based studies take wake production into account (Shaw & Patton Reference Shaw and Patton2003). Concerning the constants (
$\beta _p^{LES}$ and 
$\beta _d^{LES}$) and their values, the uncertainty involved in LES-based modelling is much lower, since only one constant and two values (8/3 or 2) are used, and in the case of the subgrid-scale dynamic model (e.g. Yue et al. Reference Yue, Meneveau, Parlange, Zhu, Kang and Katz2008; Silva Lopes et al. Reference Silva Lopes, Palma and Viana Lopes2013, Reference Silva Lopes, Palma and Piomelli2015; Yan et al. Reference Yan, Huang, Miao, Cui and Zhang2017) it does not even apply, because there is no equation for TKE of the subgrid scales (
$e$). Concerning the uncertainty, even arbitrariness, as regards the constants in LES-based modelling, it is much lower, since only one constant (
$\beta _d^{LES}$) and two values (8/3 or 2) are used, or, as is the case of the subgrid-scale dynamic model (e.g. Yue et al. Reference Yue, Meneveau, Parlange, Zhu, Kang and Katz2008; Silva Lopes et al. Reference Silva Lopes, Palma and Viana Lopes2013, Reference Silva Lopes, Palma and Piomelli2015; Yan et al. Reference Yan, Huang, Miao, Cui and Zhang2017), this does not even apply, because there is no equation for TKE of the sub-grid scales (
$e$).
Our goal in this work was to derive the most accurate equations for the effect a canopy has on a turbulent flow. The canopy-related terms in the equations, namely momentum, TKE and its dissipation rate, were approximated using a Taylor series expansion of the velocity magnitude, with a velocity scale based on the total kinetic energy. Unlike previous models (table 1), this approximation did not require any empirical coefficient. As in previous studies (Silva Lopes et al. Reference Silva Lopes, Palma and Viana Lopes2013, Reference Silva Lopes, Palma and Piomelli2015; Silva Lopes, Palma & Viana Lopes Reference Silva Lopes, Palma and Viana Lopes2016), we considered LES to be a surrogate for a well-controlled experiment and we used the time-averaged LES results to validate the decomposition and determine the highest order required for an accurate representation of the effect of the trees in each of the three equations.
There are some similarities between ours and previous studies. For instance, Belcher, Jerram & Hunt (Reference Belcher, Jerram and Hunt2003) and Belcher, Finnigan & Harman (Reference Belcher, Finnigan and Harman2008) also made use of some form of decomposition in their linear versions of flow over forested regions. Our work is more similar to that of Lien et al. (Reference Lien, Yee and Wilson2005), which addressed urban canopy flows. Both our study and that of Lien et al. (Reference Lien, Yee and Wilson2005) decompose the canopy-related terms in a Taylor series and assume that the turbulent fluctuations are much smaller than the magnitude of the mean-flow velocity (weak turbulence regime). However, in this study, the velocity scale is different and validation, with a larger number of more demanding test cases, shows that the new approach does not suffer from the limitations raised by Lien et al. (Reference Lien, Yee and Wilson2005) and that it can deal with low-velocity regions and flow-developing regions.
1.1. Paper outline
Section 2 addresses the canopy-related terms in the transport equations for momentum, Reynolds stresses, TKE and its dissipation rate, which are approximated by a Taylor series expansion in § 3. The quality of the expansion and the lowest order required for an accurate evaluation of the canopy effects are analysed in § 4, on the basis of LES results for the flow over a horizontally homogeneous forest, a forest edge and a forested hill. Section 5 presents a summary of our conclusions.
2. Mathematical model
A neutrally stratified flow within a forest is governed by equations expressing the conservation of mass and momentum:
\begin{gather} \frac{\partial \langle u_i \rangle}{\partial x_i} = 0, \end{gather}
\begin{gather}\frac{\partial \langle u_i \rangle}{\partial t} + \frac{\partial (\langle u_j \rangle \langle u_i \rangle)}{\partial x_j} = - \frac{1}{\rho} \frac{\partial \langle p \rangle}{\partial x_i}+ \nu \frac{\partial^2 \langle u_i \rangle}{\partial x_j \partial x_j} - \frac{\partial \tau_{ij}}{\partial x_j} + f_{i}, \end{gather}
where 
$u_i$ is the component of the velocity in the Cartesian direction 
$x_i$, 
$p$ is the pressure, 
$\rho$ and 
$\nu$ are the standard air density and kinematic viscosity and 
$\langle \ \rangle$ denotes a spatial filtering operation. Parameters 
$\tau _{ij} = \langle u_i u_j \rangle - \langle u_i \rangle \langle u_j \rangle$ are the LES subgrid stresses. The canopy drag force per unit mass 
$f_i$ arises because filtering and derivation do not commute in the space occupied by the canopy elements. It is usually determined by
\begin{equation} f_i = - C_d \, a(z) |\langle\boldsymbol{u}\rangle| \langle u_i \rangle , \end{equation}
where 
$C_d$ is a drag coefficient (assumed constant, since form drag is much larger than viscous drag) and 
$|\langle \boldsymbol {u}\rangle | = ( \langle u_i \rangle \langle u_i \rangle )^{1/2}$ is the scalar velocity. The spatial filtering over foliage elements used to obtain (2.1) and (2.2) was formalized by Wilson & Shaw (Reference Wilson and Shaw1977) and Raupach & Shaw (Reference Raupach and Shaw1982), and reviewed in detail by Finnigan & Shaw (Reference Finnigan and Shaw2008).
After time averaging – denoted by an overbar – of the momentum equation (2.2),
\begin{equation} \frac{ \partial \overline{\langle u_i \rangle}}{\partial t} + \frac{\partial ( \overline{\langle u_j \rangle \langle u_i \rangle )}}{\partial x_j} = \frac{\textrm{D} \overline{\langle u_i \rangle}}{\textrm{D} t} = \cdots+ \mathcal{F}_{u_i}, \end{equation}the average canopy drag is found to be
\begin{equation} \mathcal{F}_{u_i} = \overline{f_{i}} = - C_d a(z) \overline{|\langle\boldsymbol{u}\rangle| \langle u_i \rangle} . \end{equation}
Henceforth, we will drop the angle brackets that denote the filtering operation when referring to time averages of filtered quantities, i.e. 
$\bar {u}_i \equiv \overline {\langle u_i \rangle }$ and 
$\overline {u'_i u'_j} \equiv \overline {\langle u_i \rangle ' \langle u_j \rangle '}$. The Reynolds stresses are 
$\overline {u'_i u'_j}$, where 
$u'_i = u_i - \bar {u}_i$ is the fluctuating part of the instantaneous velocity.
The counterpart of 
$\mathcal {F}_{u_i}$ in the transport equations for the Reynolds stresses 
$\overline {u'_i u'_j}$ is obtained in a similar manner to the other terms in these transport equations:
\begin{equation} \overline{u'_i \frac{\textrm{D} u'_j}{\textrm{D} t} + \frac{\textrm{D} u'_i}{\textrm{D} t} u'_j} = \frac{\textrm{D} \overline{u'_i u'_j}}{\textrm{D} t} = \cdots + \mathcal{F}_{u'_i u'_j} , \end{equation}where
\begin{equation} \mathcal{F}_{u'_i u'_j} = \overline{u'_i \, f_j + u'_j \, f_i} \end{equation}
is the correlation between canopy drag and velocity fluctuations. In the same way, for each component 
$\varepsilon _{ij}$ of the Reynolds stresses dissipation rate tensor,
\begin{equation} 2 \nu \overline{\frac{\partial u'_i}{\partial x_k} \frac{\textrm{D}}{\textrm{D} t} \left[ \frac{\partial u'_j}{\partial x_k} \right] + \frac{\textrm{D}}{\textrm{D} t} \left[ \frac{\partial u'_i}{\partial x_k} \right] \frac{\partial u'_j}{\partial x_k}} = \frac{\textrm{D}}{\textrm{D} t} \; \overline{2 \nu \frac{\partial u'_i}{\partial x_k} \frac{\partial u'_j}{\partial x_k}} = \frac{\textrm{D} \varepsilon_{ij}}{\textrm{D} t} = \cdots + \mathcal{F}_{\varepsilon_{ij}}, \end{equation}with
\begin{equation} \mathcal{F}_{\varepsilon_{ij}} = 2 \nu \overline{\frac{\partial u'_i}{\partial x_k} \frac{\partial f_j}{\partial x_k} + \frac{\partial u'_j}{\partial x_k} \frac{\partial f_i}{\partial x_k}}\end{equation}
being responsible for the enhancement of dissipation by the correlation between the gradient of the canopy drag and the turbulent strain rate. Finally, the terms that correspond to the action of the canopy drag on the TKE (or 
$k$) and on its dissipation rate (
$\varepsilon$) are obtained by contracting (2.7) and (2.9) for 
$\langle u'_i u'_i \rangle$ and 
$\varepsilon _{ii}$:
\begin{gather} \mathcal{F}_k = \overline{u'_i \, f_i} = - C_d \, a(z) \overline{|\boldsymbol{u}| u_i u'_i}, \end{gather}
\begin{gather}\mathcal{F}_{\varepsilon} = 2 \nu \overline{\frac{\partial u'_i}{\partial x_k} \frac{\partial f_i}{\partial x_k}} = -2 \nu C_d \overline{\frac{\partial}{\partial x_k} \left[ a(z) |\boldsymbol{u}| u_i \right] \frac{\partial u'_i}{\partial x_k}}. \end{gather}
Because 
$f_i$ (equation (2.3)) results from the spatial filtering, 
$\mathcal {F}_k$ and 
$\mathcal {F}_{\varepsilon }$ only represent the canopy action at larger scales. This means that 
$\mathcal {F}_k$ cannot account for the increased turbulent production caused by enhanced strain of the flow in the space between leaves – wake production. However, such small-scale turbulence should dissipate quickly (over a short distance) and contribute little to TKE transport (Raupach & Shaw Reference Raupach and Shaw1982; Wilson Reference Wilson1988; Shaw & Patton Reference Shaw and Patton2003).
3. Weak turbulence regime
Since the canopy-related terms (equations (2.3), (2.10) and (2.11)) depend on the velocity magnitude, Lien & Yee (Reference Lien and Yee2005) introduced an approximation for 
$|\boldsymbol {u}|$ based on perturbation theory:
\begin{equation} |\boldsymbol{u}| = U \sqrt{1 + 2 \frac{u'_{{\textit{s}}}}{U} + \frac{u'_i u'_i}{U^2}}, \end{equation}
where 
$U = ( U_i \, U_i )^{1/2}$ is the magnitude of the mean velocity (
$U_i = \bar {u}_i$ ) and 
$u'_{{\textit{s}}}$ and 
$u'_{{\textit{n}}}$, given by
\begin{equation} u'_{{\textit{s}}} \equiv u'_i \frac{U_i}{U} \quad \text{and} \quad u'_{n} \equiv \left( u'_i u'_i - u'_{{\textit{s}}} u'_{{\textit{s}}} \right)^{1/2}, \end{equation}are the velocity fluctuations along the streamwise and normal directions defined by the mean flow. Using a series expansion of the square root in (3.1), Lien & Yee (Reference Lien and Yee2005) expressed the canopy-related terms in the transport equations of a RANS model as the sum of time-averaged second-order velocity moments. The expansion about the mean flow state assumes that the turbulent fluctuations are much smaller than the magnitude of the mean-flow velocity – weak turbulence regime (WTR):
\begin{equation} u'_i u'_i \ll U^2. \end{equation}This approach can be applied near the surface, because both the mean velocity and the fluctuations approach zero. However, in a region of separated flow, as on the lee side of a hill covered with vegetation, the turbulent fluctuations can be large and the mean velocity small, or even null, which breaks the WTR assumption (3.3).
Here, we propose to redefine the velocity scale (3.1), expanding 
$|\boldsymbol {u}|$ around a velocity 
$Q$, derived from the total kinetic energy (mean flow + turbulence):
\begin{equation} |\boldsymbol{u}| = Q \sqrt{1 + 2 \frac{U u'_{{\textit{s}}}}{Q^2} + \frac{u'_i u'_i - 2 k}{Q^2}}, \end{equation}
with 
$Q = ( U^2 + 2k )^{1/2}$. We define 
$\xi$ and 
$\eta$ such that
\begin{equation} \xi = \frac{U u'_{{\textit{s}}}}{Q^2} \quad \text{and} \quad \eta = \frac{u'_i u'_i - 2 k}{Q^2} , \end{equation}which will be used as the small parameters required to perform the Taylor series expansion of the square root:
\begin{equation} \sqrt{1 + 2 \xi + \eta} \approx 1 + \xi + \frac{\eta - \xi^2}{2} + \frac{\xi^3 - \xi \eta}{2} + \frac{6 \xi^2 \eta - \eta^2 - 5 \xi^4}{8} + \cdots . \end{equation} In the case of Lien & Yee (Reference Lien and Yee2005) (3.1), the expansion of the fluctuating term has a non-zero average, yielding an expansion around zero, fluctuations and a poor quality of series convergence. In (3.4), we reorganized the argument of the square root to have a fluctuating term around zero and a greater velocity scale (
$Q > U$). Expanding with the rescaled fluctuating term leads to a better convergence of the Taylor series.
An indication of the validity of the approximation can be inferred from the joint probability distribution of velocity fluctuations in a LES of a horizontally homogeneous canopy. In a statistically steady-state regime, the velocity fluctuations are well described by a probability distribution that is a function only of the distance to the ground 
$z$. The probability distribution is given by
\begin{equation} P (\boldsymbol{u}|\boldsymbol{r}) = \lim_{T \to \infty}\frac{1}{T} \int_t^{t + T} \delta ( \boldsymbol{u'} (\boldsymbol{r},t') - \boldsymbol{u} ) \, \mathrm{d} t' \end{equation}
and any average of the fluctuations at a point in space, 
$\boldsymbol {r}$, can be obtained by a phase-space average:
\begin{equation} \bar{A} (\boldsymbol{r}) = \int A (\boldsymbol{u'}) P (\boldsymbol{u'}|\boldsymbol{r}) \,\mathrm{d} \boldsymbol{u'}, \end{equation}
where 
$A$ is a generic function of the velocity.
A Taylor series expansion of the operator 
$A (\boldsymbol {u'})$ is valid if the higher orders of the series expansion are a good approximation in the regions where the probability weight is relevant. Figure 1 shows the joint probability distribution of the streamwise and normal components of the velocity fluctuations (obtained from LES results, as in § 4.2) at the location where the turbulence is strongest. We can observe that the probability weight is concentrated near the origin, where the instantaneous velocity is equal to the mean velocity. This measure suggests the approximation based on small velocity fluctuations (WTR) is valid and shows that the distribution is concentrated in a region inside the average fluctuation 
$(u'_{{\textit{s}}}, u'_{{\textit{n}}} )$ and that the reduced variables 
$\xi$ and 
$\eta$ (equation (3.5)) will be small even when the mean velocity 
$U$ is small.
Figure 1. Joint probability density function of the streamwise and normal velocity components 
$( u'_{{\textit{s}}} , u'_{{\textit{n}}})$ of the flow inside a horizontally homogeneous forest.
3.1. Transport of momentum
Although (2.3) has been used to model the canopy drag with both the LES (Shaw & Schumann Reference Shaw and Schumann1992; Kanda & Hino Reference Kanda and Hino1994; Dwyer et al. Reference Dwyer, Patton and Shaw1997; Patton et al. Reference Patton, Shaw, Judd and Raupach1998) and RANS (Wilson Reference Wilson1988; Svensson & Häggkvist Reference Svensson and Häggkvist1990; Green Reference Green1992) methods, it does not yield the same time-averaged force since 
$|\boldsymbol {\bar {u}}| \bar {u}_i \leq \overline {|\boldsymbol {u}| u_i}$, as signalled by Su et al. (Reference Su, Shaw, Paw U, Moeng and Sullivan1998) and Finnigan & Shaw (Reference Finnigan and Shaw2008), who suggested using a lower 
$C_d$ with LES.
Here, we use the WTR assumption and the expansion (3.6) to propose an alternative formulation for (1.1) in the RANS framework, in order to match the LES results. The canopy drag (2.5) can be expressed as a sum of terms with a decreasing order of magnitude:
\begin{equation} \mathcal{F}_{U_{i}}=\lim_{n\to+\infty}S_{U_{i}}(n), \end{equation}
where the 
$n$-truncated Taylor expansion,
\begin{equation} S_{U_i} (n) = - C_d \, a(z) \sum_{j=0}^n s_{U_i} (\,j), \end{equation}is generated by the expansion of the square root (3.6):
\begin{equation} \frac{\overline{|\boldsymbol{u}|\boldsymbol{u}}}{Q}=U_{i}\overline{\sqrt{1+2\xi+\eta}}+\overline{u_{i}\sqrt{1+2\xi+\eta}} . \end{equation}Collecting order by order, the first four terms of (3.10) are
\begin{gather} s_{U_i} (0) = Q U_i, \end{gather}
\begin{gather}s_{U_i} (1) = 0, \end{gather}
\begin{gather}s_{U_i} (2) = \frac{U}{Q} \left( \overline{u'_{{\textit{s}}} u'_i} - \frac{U \, U_i}{2 Q^2} \, \overline{u'_{{\textit{s}}} u'_{{\textit{s}}}} \right) = \frac{U_j}{Q} \left( \delta_{ik} - \frac{U_i \, U_k}{2 Q^2} \right) \overline{u'_j u'_k}, \end{gather}
\begin{align}
s_{U_i} (3) & = \frac{1}{2 Q} \left( \overline{u'_i u'_j u'_j} - \frac{U^2}{Q^2} \overline{u'_{{\textit{s}}} u'_{{\textit{s}}} u'_i} \right)\nonumber\\
& = \frac{1}{2 Q} \left( \delta_{ij} - \frac{U_i \, U_j}{Q^2} \right) \left( \delta_{kl} - \frac{U_k U_l}{Q^2} \right) \overline{u'_j u'_k u'_l}. \end{align} The zeroth-order term (3.12a) agrees with the common RANS canopy drag (1.1) when the TKE is much lower than the kinetic energy of the mean field (
$k \ll U^2$, 
$Q \approx U$). However, there is numerical (see figure 1 of Silva Lopes et al. (Reference Silva Lopes, Palma and Viana Lopes2013)) and experimental (see figure 3 of Cescatti & Marcolla (Reference Cescatti and Marcolla2004)) evidence that the ratio between turbulent (
$k$) and mean field kinetic (
$1/2 U^2$) energies could reach between 10 and 
$100\,\%$. The differences emerge when higher-order terms are considered: for instance, the second-order term (3.12c) of the WTR expansion predicts an average finite drag normal to the flow direction if the two velocity components are correlated, i.e. the Reynolds shear stress 
$\overline {u'_{{\textit{s}}} u'_{{\textit{n}}}}$ is not null. This cannot be predicted by (1.1), the basis of the RANS canopy models. Although the normal component is one order of magnitude smaller than the streamwise component (as shown by the LES results in § 4.2), this is a major conceptual difference between this formulation and the common RANS approach (1.1).
3.2. Transport of TKE
Following the same procedure as in § 3.1, a series expansion can also be obtained for the rate of work done by the canopy drag against the turbulent fluctuations (2.10):
\begin{equation} S_k (n) = - C_d \, a(z) \sum_{j=0}^n s_k (\,j), \end{equation}where
\begin{gather} s_k (0) = 2 Q k , \end{gather}
\begin{gather}s_k (1) = 0, \end{gather}
\begin{gather}s_k (2) = \frac{U^2 (Q^2 - k)}{Q^3} \overline{u_{{\textit{s}}}'^2}, \end{gather}
\begin{gather}s_k (3) = \frac{U (3 Q^2 - 2 k)}{2 Q^3} \overline{u'_{{\textit{s}}} u'_i u'_i} - \frac{U^5}{2 Q^5} \overline{u_{{\textit{s}}}'^3}, \end{gather}
\begin{gather}s_k (4) = \frac{U^4 ( 2 Q^2 - 5 k )}{4 Q^7} \overline{u_{{\textit{s}}}'^4} - \frac{U^2 (2 Q^2 - 3 k)}{2 Q^5} \overline{u_{{\textit{s}}}'^2 ( u'_i u'_i - 2 k)}+ \frac{2 Q^2 - k}{4 Q^3} \left[ \overline{(u'_i u'_i)^2} - 4 k^2 \right]. \end{gather} If the TKE is much lower than the kinetic energy of the mean field (
$k \ll U^2$, 
$Q \approx U$) and isotropy is assumed, the second-order approximation is identical to the RANS 
$k\text {--}\varepsilon$ canopy models (1.2) without wake production (
$\beta _p = 0$). This is a consequence of the assumptions of our mathematical model (§ 2) and agrees with the LES results obtained by Silva Lopes et al. (Reference Silva Lopes, Palma and Viana Lopes2013), in whose work the canopy was also modelled as a homogeneous porous medium. Wake production is due to horizontal heterogeneity of the vegetation and exists at scales much smaller than the grid sizes currently used in the simulation of atmospheric flows. This is different from canopy-resolving LES as applied to rod or urban canopy flows (e.g. Novak et al. Reference Novak, Warland, Orchansky, Ketler and Green2000; Castro, Cheng & Reynolds Reference Castro, Cheng and Reynolds2006; Cava & Katul Reference Cava and Katul2008; Bou-Zeid et al. Reference Bou-Zeid, Overney, Rogers and Parlange2009; Patton et al. Reference Patton, Horst, Sullivan, Lenschow, Oncley, Brown, Burns, Guenther, Held and Karl2011; Poggi, Katul & Vidakovic Reference Poggi, Katul and Vidakovic2011; Blackman et al. Reference Blackman, Perret, Calmet and Rivet2017). Truncating the series to second order requires only quantities that can be obtained from RANS models. Third-order terms can also be included if an approximation for the third-order correlations is used.
3.3. Transport of the dissipation rate of TKE
Finally, an approximation for the effect of the canopy drag on the transport of TKE dissipation rate (2.11) was also derived:
\begin{equation} S_{\varepsilon} (n) = - \nu \, C_d \sum_{j=0}^n s_{\varepsilon} (\,j) . \end{equation}
The expansion of the product between the gradient of the canopy drag and the turbulent strain rate is more complex than that for momentum or TKE transport. For the sake of convenience (§ 4.4), the series terms were split into large-scale contribution, 
$S_{\varepsilon }^{{\textit{L}}} (n)$, which depends on velocity correlations, and small-scale contribution, 
$S_{\varepsilon }^{{\textit{s}}} (n)$, which depends on statistics that include the turbulent strain rate:
\begin{equation} S_{\varepsilon} (n) = S_{\varepsilon}^{{\textit{L}}} (n) + S_{\varepsilon}^{{\textit{s}}} (n) = - \nu \, C_d \left( \sum_{j=0}^n s_{\varepsilon}^{{\textit{L}}} (\,j) + \sum_{j=0}^n s_{\varepsilon}^{{\textit{s}}} (\,j) \right) ,\end{equation}where
\begin{gather} s_{\varepsilon}^{{\textit{L}}} (0) = s_{\varepsilon}^{{\textit{L}}} (1) = 0, \end{gather}
\begin{gather}s_{\varepsilon}^{{\textit{L}}} (2) = \frac{\partial}{\partial x_k} \left[ a(z) Q \left( \frac{U_i \, U_j}{Q^2} + \delta_{ij} \right) \right] \frac{\partial \overline{u'_i u'_j}}{\partial x_k} , \end{gather}
\begin{gather}s_{\varepsilon}^{{\textit{L}}} (3) = \frac{\partial}{\partial x_l} \left[ a(z) \, \frac{U_i}{Q} \left( \delta_{jk} - \frac{U_j \, U_k}{3 Q^2} \right) \right] \frac{\partial \overline{u'_i u'_j u'_k}}{\partial x_l} \end{gather}and
\begin{gather} s_{\varepsilon}^{{\textit{s}}} (0) = s_{\varepsilon}^{{\textit{s}}} (1) = 0, \end{gather}
\begin{gather}s_{\varepsilon}^{{\textit{s}}} (2) = 2 a(z) \, Q \left( \delta_{ij} + \frac{U_i \, U_j}{Q^2} \right) \overline{\frac{\partial u'_i}{\partial x_k} \frac{\partial u'_j}{\partial x_k}}, \end{gather}
\begin{gather}s_{\varepsilon}^{{\textit{s}}} (3) = 2 a(z) \left( \frac{U_i \delta_{jk} + U_j \delta_{ik} + U_k \delta_{ij}}{Q} - \frac{U_i U_j U_k}{Q^3} \right) \overline{u'_i \frac{\partial u'_j}{\partial x_l} \frac{\partial u'_k}{\partial x_l}} . \end{gather} As with TKE transport, if the TKE is much lower than the kinetic energy of the mean field (
$k \ll U^2$, 
$Q \approx U$) and isotropy of both the Reynolds stresses and TKE dissipation rate is assumed, the second-order small-scale approximation becomes similar to conventional RANS 
$k\text {--}\varepsilon$ canopy models without wake production (equation (1.3)). As will be shown in § 4.4, the magnitude of the small-scale part overshadows the large-scale part, since dissipation is mainly effected at the smaller scales. Also, if isotropy is assumed, the second-order approximation uses quantities available in RANS models. The correlations required by the third-order small-scale term are difficult to approximate.
4. Numerical experiments and validation
Large-eddy simulations were performed (table 2) to assess the accuracy of the proposed model (equations (3.12), (3.14), (3.17) and (3.18)), that is, the function of the order of approximation in the description of canopy effects on the transport of momentum, TKE and its dissipation rate. Three different flows – over a horizontally homogeneous canopy, a forest edge and a forested hill – were considered on the basis of studies by Shaw & Schumann (Reference Shaw and Schumann1992), Yang et al. (Reference Yang, Morse, Shaw and Paw U2006) and Dupont, Brunet & Finnigan (Reference Dupont, Brunet and Finnigan2008), respectively. These flows have already been considered by Silva Lopes et al. (Reference Silva Lopes, Palma and Viana Lopes2013), because they test the accuracy of the model for flows from the simplest to one with pressure gradients, streamline curvature and separation. The results obtained by Silva Lopes et al. (Reference Silva Lopes, Palma and Viana Lopes2013) and their set of optimized coefficients, identified as reference RANS 
$k\text {--}\varepsilon$, will serve as an additional source for comparing the performance of the decomposition.
Table 2. Numerical parameters. Here 
$U_b$, bulk velocity; 
$LAI = \int _0^h a(z)\, \mathrm {d} z$, leaf area index; 
$f_{z_{\textit{max} }}$, maximum expansion factor of the grid in the vertical direction.
4.1. Numerical model
The numerical model is composed of (2.1) and (2.2) and the subgrid stresses, modelled by an eddy-viscosity assumption:
\begin{equation} \tau_{ij} - \tfrac{1}{3} \delta_{ij} \tau_{kk} = -2 \nu_{{\boldsymbol{\mathsf{t}}}} S_{ij} = -2 \left( C \Delta \right)^2 | S | S_{ij}, \end{equation}
where 
$\varDelta = ( \Delta x \, \Delta y \, \Delta z )^{1/3}$ is the filter size, 
$S_{ij} = (\partial u_i / \partial x_j + \partial u_j / \partial x_i ) / 2$ is the resolved strain-rate tensor and 
$| S | = (2 S_{ij} S_{ij})^{1/2}$ is its magnitude. The coefficient 
$C$ was obtained using a combination of the Meneveau, Lund & Cabot (Reference Meneveau, Lund and Cabot1996) Lagrangian dynamic model with the wall-damped Mason & Thomson (Reference Mason and Thomson1992) Smagorinsky model near the surface. The two models merged below 
$0.75 h$ and 
$C_{{{\textit {S}}}} = 0.16$ was used in the Smagorinsky model. Further details of the implementation can be found in Silva Lopes, Palma & Castro (Reference Silva Lopes, Palma and Castro2007) and Silva Lopes et al. (Reference Silva Lopes, Palma and Viana Lopes2013). Periodic conditions were used in the streamwise and spanwise directions. The surface stress was determined using the Marusic, Kunkel & Porté-Agel (Reference Marusic, Kunkel and Porté-Agel2001) model, while the top of the domain used a free-slip condition. At each time step, a uniform (in space) streamwise pressure gradient was determined to balance the canopy drag, resulting in a constant flow rate (momentum balance).
Numerical uncertainty has been assessed previously (Silva Lopes et al. Reference Silva Lopes, Palma and Viana Lopes2013), on the basis of the streamwise velocity normalized with the volume-averaged streamwise velocity component (
$U_b$), and the TKE and Reynolds shear stress (
$u'_{{{\textit{s}}}}u'_{{{\textit{n}}}}$). Both of these were normalized with the friction velocity (
$u_*$), calculated using the total shear stress (resolved and subgrid) at the canopy top. Calculations in three numerical grids of increasing resolution, of which the finest was identical to the one in this study, showed that 
$U/U_b$ and 
$u'_{{\textit{s}}}u'_{{\textit{n}}}/u^2_*$ were independent of grid resolution. In contrast, 
$k/u^2_*$ above the canopy increased nearly 
$10\,\%$ in peak value between the coarsest grid (
$96\times 48\times 48$) and finest grid (
$192\times 96\times 98$) results. The contribution of the subgrid stress to the peak Reynolds shear stress, at treetop level, decreased from 
$4\,\%$ with the coarsest grid to 
$2\,\%$ with the finest. The grid refinement study was performed for the horizontally homogeneous canopy, while in the other cases, the resolution was limited by computational resource availability. It was, nevertheless, considered appropriate for this study.
4.2. Transport of momentum
Figure 2 shows the WTR approximation of the canopy drag along the streamwise and normal directions in the horizontally homogeneous forest, where there is no streamline curvature, and a comparison with LES results. Regardless of the direction, the magnitude of the series terms decreases with the order and the sum approaches the LES results (figure 2a,b), suggesting the convergence of the series. The WTR expansions correctly predicted both drag components, with the largest terms being based on the streamwise velocity and on the Reynolds shear stress 
$\overline {u' w'}$. The zeroth order was sufficient to improve the standard RANS approximation, which underpredicted streamwise drag by 
$25\,\%$ (figure 2c).
Figure 2. Approximations of the canopy drag along the streamwise and normal directions (
$s$ and 
$n$) in the horizontally homogeneous forest. 
$(a)$ The WTR series terms up to the third order (
$s_{U_i}$ (3.12)) and LES results (
$\mathcal {F}_{U_i}$ (2.5)); 
$(b)$ WTR summations up to the third order (
$S_{U_i}$ (3.10)) and LES results; 
$(c)$ reference RANS 
$k\text {--}\varepsilon$ (
$\mathcal {F}^{k-\varepsilon }_{U_i}$ (1.1)) and LES results (
$\mathcal {F}_{U_i}$ (2.5)). Note that horizontal scales are magnified 
$20\times$ in the case of quantities referring to the normal direction, 
$n$.
In the case of the forest edge and the forested hill, with streamline curvature, it is useful to analyse the drag force in terms of its streamwise and normal components:
\begin{equation} \mathcal{F}_{U_s} (n)= S_{U_i} (n) \frac{U_i}{U} \quad \text{and} \quad \mathcal{F}_{U_n} (n)= \left[ \left| \boldsymbol{S_U} (n) \right|^2 - \mathcal{F}_{U_s}^2(n) \right]^{1/2} , \end{equation}
and to study canopy force modelling by verifying the correlation of the WTR expansion up to the third order (3.10) with the canopy drag determined in the LES simulations at each point within the canopy layer (figures 3 and 4). In the case of the streamwise component, a perfect match is already achieved by the second-order WTR approach (figure 3b,c), whereas the standard RANS presented deviations of up to 
$20\,\%$ (figure 3d). The normal component was around 
$20{\times }$ smaller than the streamwise one (as in the horizontally homogeneous canopy flow) and the convergence was slower, especially in the case of the forested hill (figure 4). However, the dispersion observed in the correlation was, in either case, restricted to specific regions of the physical domain (the leading edge of the forest and the lee side of the forested hill, inside the separated flow region).
Figure 3. Correlations of approximations of the streamwise component of the canopy drag in the forest edge and forested hill flows. (a–c) The WTR summations up to the third order (
$S_{U_s}$ (3.10)) and LES results (
$\mathcal {F}_{U_s}$ (2.5)); (d) reference RANS 
$k\text {--}\varepsilon$ (
$\mathcal {F}^{k-\epsilon }_{U_s}$ (1.1)) and LES results.
In summary and on the basis of the three flow cases, we conclude that accurate modelling of the canopy drag based on the WTR expansion required a second-order approximation of the streamwise component and a fourth-order approximation of the normal component. This may be compared with the standard RANS approach, which was less accurate than the zeroth-order WTR for the streamwise component and does not predict any drag in the normal direction.
4.3. Transport of TKE
As with canopy drag, there is evidence of convergence of the WTR approximation of the effect on TKE transport in the horizontally homogeneous canopy flow, since the magnitude of the terms decreased with order and their sum approached the LES results (figure 5a,b). However, convergence was slower than for canopy drag (figure 2) and the third order was required to capture 
$95.8\,\%$ of the LES results, compared to 
$83.5\,\%$ in the case of the second order (figure 5b). Also, the WTR approximation could not improve on the reference RANS model prediction of Silva Lopes et al. (Reference Silva Lopes, Palma and Viana Lopes2013):
\begin{equation} \mathcal{F}_k^{k-\varepsilon} = -C_d \, a(z) \, \beta_d |\boldsymbol{U}| k, \end{equation}
in which the constant 
$\beta _d = 4.0$ was calibrated with this flow (figure 5c).
Figure 5. Approximations of the contribution of canopy drag to TKE transport in the horizontally homogeneous forest. (a) The WTR series terms up to the fourth order (
$s_k$ (3.14)) and LES results (
$\mathcal {F}_k$ (2.10)); (b) WTR summations up to the fourth order (
$S_k$ (3.13)) and LES results; (c) reference RANS 
$k\text {--}\varepsilon$ (
$\mathcal {F}_k^{k-\varepsilon }$ (4.3)) and LES results.
In the more complex cases, the flows over the forest edge and the forested hill, the findings of the horizontally homogeneous canopy flow were confirmed, with the third order providing a good approximation and the fourth order matching the LES results (figure 6c,d). In these cases, the WTR approximation improved the reference RANS model predictions (figure 6e), which failed mainly at the edge of the forest and in the separation on the lee side of the hill. However, we note that the third-order term requires velocity correlations of the same order.
Figure 6. Correlations of approximations of the contribution of canopy drag to TKE transport in the forest edge and forested hill flows. (a–d) The WTR summations up to the fourth order (
$S_k$ (3.13)) and LES results (
$\mathcal {F}_k$ (2.10)); (e) reference RANS 
$k\text {--}\varepsilon$ (
$\mathcal {F}_k^{k-\varepsilon }$ (4.3)) and LES results (
$\mathcal {F}_k$ (2.10)).
4.4. Transport of dissipation rate of TKE
The WTR approximation of the effect of canopy drag on the transport of the dissipation rate of TKE includes only second- and higher-order terms (equations (3.17) and (3.18)). As in the transport of momentum or TKE, the magnitude of the third-order term was smaller than the second-order term and their sum approached the LES result in the horizontally homogeneous canopy flow (figure 7a,b), indicating convergence of the series. The second-order term of the series, (3.17b) and (3.18b), the first non-zero terms of the series, yielded a result within 
$83.8\,\%$ of the peak value (figure 7b). The LES results were mimicked by either the WTR third-order approximation (
$99.6\,\%$) or the reference RANS model prediction of Silva Lopes et al. (Reference Silva Lopes, Palma and Viana Lopes2013):
\begin{equation} \mathcal{F}_{\varepsilon}^{k-\varepsilon} = -C_d a(z) C_{\varepsilon 5} \beta_d |\boldsymbol{U}| \varepsilon, \end{equation}
in which the constant 
$C_{\varepsilon 5} = 0.9$, such as 
$\beta _d$ (4.3), was calibrated also with this flow (figure 7c).
Figure 7. Approximations of the contribution of canopy drag to the transport of the TKE dissipation rate in the horizontally homogeneous forest. (a) The WTR series terms up to the third order (
$s_{\varepsilon }$ (3.17) and (3.18)) and LES results (
$\mathcal {F}_{\varepsilon }$ (2.11)); (b) WTR summations up to the third order (
$S_{\varepsilon }$ (3.16)) and LES results (
$\mathcal {F}_{\varepsilon }$ (2.9)); (c) reference RANS 
$k\text {--}\varepsilon$ (
$\mathcal {F}_{\varepsilon }^{k-\varepsilon }k$ (1.3)) and LES results (
$\mathcal {F}_{\varepsilon }$ (2.9)).
For the forest edge and forested hill flows, the second-order WTR approximation was sufficient to improve the correlation of the reference RANS model, which showed considerable dispersion (figure 8). However, we note that the second-order WTR approximation requires knowledge of turbulent strain rate correlations and the third-order approximation requires third-order statistics including the turbulent strain rate and velocity fluctuations.
Figure 8. Correlations of approximations of the contribution of canopy drag to TKE dissipation rate transport in the forest edge and forested hill flows. (a) Second-order and (b) third-order WTR summations (
$S_{\varepsilon }$ (3.16)) and LES results (
$\mathcal {F}_{\varepsilon }$ (2.11)); (c) reference RANS 
$k\text {--}\varepsilon$ (
$\mathcal {F}_{\varepsilon }^{k\text {--}\varepsilon }$ (1.3)).
Additionally, in assessing the accuracy of the WTR approximation to the effect of canopy drag on the transport of TKE dissipation rate, we must bear in mind that this involves mainly small scales, which can be shown by the dependence of the WTR approximation on the turbulent strain rate correlations (3.18c). This is the opposite of the effect on the transport of momentum and TKE, which depends on velocity correlations and concerns mainly large scales. As such, a grid refinement study was performed on the horizontally homogeneous canopy flow, to check whether the accuracy of the WTR approximation varied with resolved dissipation. Three different grid resolutions were considered: low resolution (LR; 
$136 \times 68 \times 70$), medium resolution (MR; 
$192 \times 96 \times 98$) and high resolution (HR; 
$272 \times 136 \times 138$ nodes); i.e. the grid space was nearly two times larger in the LR grid than in the HR grid. As expected, the small-scale part of the WTR approximation (
$s_{\varepsilon }^{{\textit{s}}}$) increased with the grid resolution (figure 9a), whereas the large-scale part (
$s_{\varepsilon }^{{\textit{L}}}$), already properly resolved on the LR grid, remained almost constant (figure 9b). Nevertheless, there was always a good agreement between the WTR approximation determined using resolved LES fields and the LES results on the same grid.
Figure 9. WTR approximations of the contribution of canopy drag to TKE dissipation rate transport in the horizontally homogeneous forest with different grid resolutions (LR, MR and HR). (a) Second- and third-order small-scale series terms (
$s_{\varepsilon }^{{\textit{s}}}$ (3.18)); (b) Second- and third-order large-scale series terms (
$s_{\varepsilon }^{{\textit{L}}}$ (3.17)). Note that the horizontal scale is magnified 6
$\times$ in the case of the third-order terms.
As regards the grid refinement study (figure 9), if a major part of the dissipation were to be resolved, we would expect that: (a) the agreement between the WTR approximation and LES results would persist and (b) the small-scale part of the WTR approximation would overshadow the large-scale part (cf. figure 9a with figure 9b). This latter consideration can be used in the construction of a RANS model, based solely on the small-scale part. We note, however, that if the grid resolution of the LES was increased to approach the dissipative scales of the flow, the field-scale approach, used to model the forest canopy, should be replaced by the plant-scale approach (Yue et al. Reference Yue, Parlange, Meneveau, Zhu, van Hout and Katz2007). In this approach, the local canopy structure is taken into account instead, rather than a horizontally homogeneous drag field.
In the cases of the forest edge and forested hill flows (figure 10), the large- and small-scale contributions, determined using the resolved LES fields, were of a similar magnitude and only their sum correlated well with LES results. This occurred because the spatial resolution was not as good as it was for the horizontally homogeneous canopy flow. This is why fewer small scales and, consequently, dissipation were resolved.
Figure 10. Forest edge and forested hill. Correlation between the contribution of canopy drag to TKE dissipation rate transport predicted by the WTR expansion up to the third order (
$S_{\varepsilon }(3)$ (3.15)) and LES results (
$\mathcal {F}_{\varepsilon }$ (2.11)). Contribution of small-scale (
$S^{s}_{\varepsilon }(3)$ (3.16)) and large-scale (
$S^{L}_{\varepsilon }(3)$ (3.16)) components, and sum of both of these.
4.5. Comparison with conventional RANS models
The WTR approximations are similar to conventional RANS 
$k\text {--}\varepsilon$ canopy models if certain assumptions are made: (a) the TKE is much smaller than the kinetic energy of the mean field: 
$k \ll U^2$, 
$Q \approx U$); (b) the Reynolds stress and its dissipation rate tensors are isotropic: 
$\overline {u'_i u'_j} = 2 / 3 k \delta _{ij}$ and 
$\varepsilon _{ij} = 2 / 3 \varepsilon \delta _{ij}$; and (c) the large-scale contribution to the effect on the transport of TKE dissipation rate is much smaller than the small-scale part (see above): 
$S_{\varepsilon }^{{\textit{L}}} \ll S_{\varepsilon }^{{\textit{s}}}$, 
$S_{\varepsilon } \approx S_{\varepsilon }^{{\textit{s}}}$. Given these assumptions, the zeroth-order WTR approximation for canopy drag and the second-order approximation for the effect on the transport of TKE and its dissipation rate become
\begin{gather} S_U (0) = - C_d a(z) | \boldsymbol{U} | U_i, \end{gather}
\begin{gather}S_k (2) = C_d a(z) \left( - \tfrac{8}{3} | \boldsymbol{U} | k \right), \end{gather}
\begin{gather}S_{\varepsilon} (2) = C_d a(z) \left( - \tfrac{8}{3} | \boldsymbol{U} | \varepsilon \right); \end{gather}
i.e. 
$\beta _d = 2.67$ and 
$C_{\varepsilon 5} =1$.
The zeroth-order approximation for canopy drag (4.5) agrees with the conventional approach (1.1). However, the second-order WTR approximation appears to underestimate the effect on TKE and its dissipation rate when compared with the Silva Lopes et al. (Reference Silva Lopes, Palma and Viana Lopes2013) reference RANS model (figures 5c and 7c), in which the coefficients (
$\beta _d = 4.0$ and 
$C_{\varepsilon 5}=0.9$) were calibrated using LES. This is due, in large part, to the effect of the canopy being split into various orders in the WTR approximation, whereas conventional models try to account for it all with a single term, yielding larger coefficients (table 1). Also, the isotropy assumption, although useful in simplifying the equations, is not appropriate inside the canopy. It is sufficient, for instance, to consider the anisotropy commonly found in the Reynolds stresses of homogeneous shear flows, 
$\overline {u' u'} \approx k$, to obtain 
$\beta _d = 3.0$ instead of 
$2.67$. In summary, we find that the differences are justified and there is no sign of lack of accuracy in any of the approaches.
5. Summary and conclusions
The canopy-related terms in the transport equations for momentum, TKE and its dissipation rate were approximated using a Taylor series expansion of the velocity magnitude, with a velocity scale based on the total (mean flow and turbulent) kinetic energy. A consistent mathematical formulation was developed in a Reynolds-averaged formulation that, in the case of two-equation (
$k$–
$\varepsilon$) turbulence models, implies no additional empirical model coefficients or any uncertainties associated with determining the same. It was shown that in homogeneous canopy flows, the standard RANS formulation is a zeroth-order approximation of momentum and a second-order approximation of both TKE and its dissipation in the WTR formulation.
The conclusions drawn from LES of three canopy flows (over a horizontally homogeneous forest, a sequence of forest stands and clearings, and a forested hill) were as follows:
(i) The limitations inherent in the weak turbulence assumptions (i.e. turbulence intensity much lower than mean velocity or low TKE compared with mean kinetic energy) were overcome by reverting to a velocity scale derived from the total kinetic energy. The model is robust and provided good results, even where these assumptions are questionable, as is the case for non-homogeneous and non-equilibrium flow regions, found near forest edges and on the lee side of forested hills. In the three cases, the perturbative description was valid, leading to a fast-converging series.
(ii) For an accurate representation of the canopy effects in a Reynolds-averaged formulation, momentum in the streamwise direction must be approximated by up to second-order terms.
(iii) The present approach also shows the existence of a drag force perpendicular to the streamwise flow direction. This force, which is about 20
$\times$ lower than the force along the streamwise direction, must be approximated by second- and third-order terms.(iv) The TKE requires third-order terms. In the case of the horizontally homogeneous canopy, a
$96\,\%$ agreement with the LES results was achieved, with the first $84\,\%$ being the contributions of the zeroth- and second-order terms.(v) The modelling of the transport of the TKE dissipation rate was validated under conditions involving only a small fraction of the dissipation being resolved. Nevertheless, a comparison of results obtained with different grid resolutions indicates that the model is also valid and a third-order approximation yields a
$100\,\%$ agreement with the LES results for the horizontally homogeneous canopy. The series development also shows that the zeroth- and first-order terms are equal to zero, and the second-order term, the lowest approximation of the TKE dissipation rate, yields an accuracy of $84\,\%$ in the case of the horizontally homogeneous canopy.
The novelty of this model is that it allows for the determination of canopy effects on the transport of momentum, TKE and its dissipation rate, using a fast-converging series expansion. This, which contrasts with the dimensional arguments used in other turbulence models, may provide the insight required to develop more accurate approaches, thus improving the accuracy of canopy flow prediction in the RANS formulation. We note that although the 
$k$–
$\varepsilon$ model was considered here, the same approach can be applied to second-order Reynolds stress transport models. Results of and assumptions inherent to our model, such as not accounting for TKE wake production and volume fraction, and the uncertainty of both drag coefficient and leaf area index are expected to be minor factors in real life. In such applications, two main difficulties may be encountered: first, our conclusions were based on a comparison of the magnitude of each term of the series with the results of a well-controlled (numerical) experiment, which may differ from laboratory or field measurements; and second, there is the issue of how to translate the higher-order terms of our series into new approximations for canopy-related terms in the RANS formulation, which will depend on the turbulence model being used. These topics are worth addressing in future studies.
Acknowledgements
Whilst engaged in this study, J.V.L. and A.S.L. were research fellows for the Ciência 2007 and Ciência 2008 programmes run by the Portuguese Foundation for Science and Technology (FCT, Fundação para a Ciência e Tecnologia). The authors would like to acknowledge the financial support of both the European Commission and the Portuguese Foundation for Science and Technology through the ERA-NET PLUS NEWA – 2014, New European Wind Atlas Joint Programme (contract no. NEWA/0001/2014). This study was partly developed under a research contract with Renewable Energy Systems Ltd and Natural Power Consultants Ltd.
Declaration of interests
The authors report no conflict of interest.
Open access
