A CFD‐Based Mixing Model for Vegetated Flows

This paper provides a computational fluid dynamics (CFD)‐based modeling framework for predicting flow field, turbulence, and mixing characteristics within vegetated environments such as ponds and wetlands. The framework has been implemented within a commercial CFD code—ANSYS Fluent 19—via a set of user‐defined functions. Following the approach outlined by King et al. (2012, https://doi.org/10.1017/jfm.2012.113), the standard k‐ε turbulence closure model has been modified to capture the energy transfer at the vegetation/clear flow shear interface and within the vegetation. The implementation assumes that vegetation is vertical, but nonorthogonal flow in the horizontal plane is accounted for. Values for the drag coefficient and the mixing coefficients are estimated based on the vegetation stem diameter and density. Following Tanino and Nepf (2008, https://doi.org/10.1017/S0022112008000505), a switch has been incorporated to account for the fact that the relevant length scale changes from stem diameter to stem spacing as stem density increases. A set of model parameters is proposed, based on a reevaluation of previously published laboratory data and theoretical analysis. Five different experimental data sets are used to demonstrate that the model is able to predict mixing within fully vegetated systems and due to both vertical and horizontal shear layers. The framework was developed to provide a practical prediction tool for engineering purposes, in particular for the estimation of residence time distributions in real partially vegetated stormwater management ponds. Its implementation here within a commercial CFD package potentially facilitates application to complex pond geometries, including patches of different types of vegetation with different bulk stem diameter and density characteristics.


Introduction
Vegetated flows are complex and of interest in both natural and engineered systems, where the presence of vegetation can significantly affect the flow field (Nepf, 2012). At the river scale, vegetation can reduce conveyance and trap pollutants (Darby, 1999). In a treatment pond, vegetation can promote short circuiting and reduce efficiency . In wetlands, vegetation is critical to ensuring treatment. Vegetation itself can absorb pollutants; it can also form a habitat for biofilms and other microorganisms that treat pollutants, and it acts to slow velocities, allowing suspended pollutants to settle (Kadlec & Wallace, 2008;Shilton, 2005). Modeling the presence of vegetation is necessary to predict velocities and other mean flow field effects (Marjoribanks et al., 2017;Nepf, 1999& Tsavdaris et al., 2013, while solute transport modeling in vegetated systems is necessary to predict mixing and treatment processes. Computational fluid dynamics (CFD) models are useful for understanding and predicting hydrodynamic processes (Tu et al., 2012). Vegetation may be taken into consideration within CFD models by either explicitly modeling each stem in the domain (e.g., Stoesser et al., 2010) or spatially averaging the effects of vegetation using a bulk scale approach. Gac (2014) developed a lattice-Boltzmann method-based approach to the stem-scale numerical simulation of rigid and flexible submerged vegetation using a large-eddy simulation (LES)-based viscosity closure. The results showed good agreement for submerged vegetation but required approximately 200 million cells for a 1.5-L volume. Okamoto and Nezu (2010) investigated flow structure and mass transport in submerged vegetated flow using stem-scale modeling and LES. Their simulation of a small array of stems reproduced mean flow, turbulence, and concentration characteristics well. Lu & Dai (2018) investigated solute transport for rigid and flexible submerged vegetation using an LES model, integrating mechanical dispersion using a random-walk model. Their rigid vegetation solute transport results showed that concentrations were underpredicted in the free water, and they concluded that canopy density has a greater impact on solute transport than vegetation rigidity. Stoesser et al. (2009) performed LES for flow through idealized submerged vegetation. They reported "fairly good" agreement with velocity and turbulence intensities and showed the formation of vortex rolls above the vegetation interface. However, this application was limited to 16 submerged circular cylinders.
The explicit modeling of stems within a real vegetated system provides two major challenges. First, it is computationally expensive, as it requires mesh cell sizes to be of the order of the scale of the vegetation (e.g., a few millimeters) compared with the domain (e.g., tens of meters), resulting in a large number of cells. None of the examples cited above has modeled more than a small fraction of a full-scale pond volume. Second, while we may know what species grow in ponds, and have some idea of their generic characteristics, a stem-scale modeling approach would require explicit knowledge of both the location of each plant element and their functional characteristics (e.g., flexure) at the time of modeling.
The application of a drag force as a momentum sink within CFD models to represent vegetation provides an alternative to the stem-scale approach (Tsavdaris et al., 2013). Marjoribanks et al. (2017) applied a mass flux scaling algorithm to model the effects of vegetation, on a cell-by-cell basis. This used a momentum sink approach and solid volume fraction to account for porosity. While providing an approach that can describe patch-scale differential advection, the mass flux scaling algorithm ignores stem-scale vegetation-induced turbulence and mixing. Yan et al. (2017) and  both employed a similar bulk-scale approach, representing the vegetation with uniform drag and then specifying a fixed magnitude for the turbulent diffusion coefficient within the vegetation to describe the velocity field.
A more elegant approach would be to quantify the magnitude of the in-vegetation turbulence and mixing from the contributing physical flow and turbulence processes. This physically based approach would include individual estimates of the four stem-scale processes that contribute to mixing within vegetation: turbulent diffusion, mechanical dispersion (Tanino & Nepf, 2008), vortex trapping, and secondary wake dispersion (White & Nepf, 2003). These processes are specifically caused by the interaction of a stem with the flow. When the stem is not explicitly modeled, these processes are not represented within the model. The aim of this paper is to develop a modeling framework that realistically models flows and solute transport in vegetated systems by representing stem-scale processes as bulk-scale phenomena. A key requirement of the modeling framework is that it should be capable of describing flow and solute transport within full-scale storm water management ponds that include patches of vegetation with different bulk stem diameter and density characteristics.

Literature Review
The following sections summarize current literature relevant to the development of the proposed bulk-scale modeling framework. Identified knowledge gaps lead to the specific objectives for the current research, which are outlined in section 2.5.

Bulk-Scale Modeling of Vegetation Using CFD
Within a CFD model, the time-averaged Navier-Stokes equations are solved for continuity, and conservation of momentum, where ρ is density, u i is instantaneous velocity averaged over the cell volume, i = {x, y, z} directions, p is pressure, τ ij the stress tensor, and F i any external forces acting on the flow (Versteeg & Malalasekera, 2007). Reynolds averaging is typically used to solve these equations in applied contexts with a turbulence closure model.
Following the standard drag force approach (e.g., Kadlec, 1990;Nepf, 1999), drag caused by the vegetation can be modeled by where ϕ is solid volume fraction, C D is bulk drag coefficient, and a is frontal-facing area (the area of stems perpendicular to flow per unit plan area). The term (1 − ϕ) −1 takes into account the volume of water excluded by the stems when calculating driving forces (King et al., 2012). Assuming a cylinder array, solid 10.1029/2018WR023628

Water Resources Research
volume fraction is typically linked to frontal-facing area and stem diameter through ϕ = adπ/4, where d is stem diameter.
k-ε-based turbulence models, derived from the standard k-ε turbulence model of Launder and Spalding (1974), are the most widely used turbulence closures. While unmodified k-ε models can be shown to be applicable to vegetated shear layers (e.g., Fischer-Antze et al., 2001), they are not suitable for predicting mixing within vegetation in the absence of a shear layer, as they lack a turbulence production term due to vegetation (Sonnenwald et al., 2016). Recent laboratory studies, employing both particle image velocimetry and laser Doppler anemometry measurements, have quantified the turbulent kinetic energy (TKE) budget within arrays of random emergent cylinders (e.g., Ricardo et al., 2014;Ricardo et al., 2016). These measurements have confirmed that the major source of TKE is vortex shedding from individual cylinders. Choi and Yang (2005) provide a review of several turbulence closure methods suitable for bulk-scale modeling of vegetation. Several authors have implemented first-order closures, typically k-ε based, for example, López and García (2001), Katul et al. (2004), Defina and Bixio (2005), King et al. (2012), and Liu et al. (2017). Other authors have proposed higher-order closures based on the algebraic stress model or Reynolds stress model, for example, Naot et al. (1996), Wilson and Shaw (1977), and Choi and Kang (2004).
Compared with first-order closures, second-or third-order closure models have the potential to provide improved representation of the flow field in certain situations, for example, in the presence of secondary circulations and where the geometry and vegetation characteristics are accurately described. However, their increased computational expense may not always be justified and may not be practical for modeling more geometrically complex systems, such as a ponds.
The King et al. (2012) k-ε-based turbulence closure was chosen as the basis for our proposed framework for modeling solute transport within vegetation. The key motivations for selecting this model were its additional transport term for wake-scale TKE, resolving weaknesses in the dissipation of TKE, and that it shares much of its theory with Tanino and Nepf (2008), making it ideally suited to adaptation for solute transport modeling in vegetated flows. It is also conceptually simple, is numerically straightforward to implement within existing CFD codes, and has been shown to reproduce turbulence characteristics for vegetation in a range of flow conditions.

The King et al. (2012) k-ε Turbulence Model
The King et al. (2012) model represents the effects of vegetation by separating TKE and turbulent dissipation into two components: shear scale and vegetation stem scale (wake scale): where k is TKE, ε is turbulent dissipation rate, subscript s indicates shear scale, and subscript w indicates wake scale. King et al. (2012) modified the standard k-ε model shear kinetic energy transport equation to include the "spectral shortcut" proposed by Finnigan (2000): where μ is the molecular viscosity, μ t is the turbulent viscosity, σ k is a standard k-ε model constant, P s is the standard k-ε model production term for TKE due to shearing, and W is the energy lost due to the spectral shortcut. W represents a direct energy cascade from shear to stem scales, that is, a large eddy hitting a stem and breaking up immediately to a smaller length scale rather than dissipating normally. King et al. (2012) proposed that the energy from W is added to the wake scale, and included a new transport equation for wake-scale TKE k w : where P w is the production of TKE due to stem wakes and W is the additional energy from the spectral shortcut.

10.1029/2018WR023628
Water Resources Research King et al. (2012), following Nepf (1999) and Tanino and Nepf (2008), assumed that the production of wakescale TKE is due to momentum dissipation, that production is equal to dissipation, and that dissipation follows the standard k-ε model formula, that is, P w = ε w and ε w ∝ k w 3/2 d −1 , giving where β p is a model constant and U is streamwise velocity. King et al. (2012) decomposed the drag coefficient into inertial and viscous components. C Dp is the inertial (or pressure) drag coefficient. Equation (7) only converts inertial losses, not viscous losses, into TKE. The method we adopt to define C Dp is described in section 3.1. King et al. (2012) defined the spectral shortcut as where β d is another model constant. The standard k-ε equation for the transport of shear-scale TKE dissipation (ε s ) is modified due to the spectral shortcut: where C ε1 , C ε2 , and σ ε are standard k-ε model constants from Launder and Spalding (1974) and C ε5 is a model constant. Wake-scale TKE dissipation is modeled directly through the following relationship: where C εD is a model constant. King et al. (2012) noted that d in equation (10) should change to s, the stemedge to stem-edge spacing, when s < d. King et al. (2012) modified the turbulent viscosity to take into account both the shear and wake scales of TKE, giving where C μ is a standard k-ε model constant from Launder and Spalding (1974) and C λ is a new model constant.
In the absence of vegetation, the King et al. (2012) model collapses to the standard k-ε model. In the absence of shear, the production of turbulence due to stems must equal its dissipation.
In total, there are five new model constants: β p , β d , C ε5 , C εD , and C λ . King et al. (2012) estimated appropriate values for these constants by fitting to their experimental data and data described in Dunn et al. (1996) to give β p = 0.2, β d = 1.0, C ε5 = 0.0, C εD = 0.28, and C λ = 0.01. Setting the value of C ε5 equal to 0.0 means that in practice the transport of shear-scale turbulent energy dissipation (equation (9)) remains unchanged from the standard k-ε model.
The original King et al. (2012) model has some limitations. Its dependency on the calibration of five constants leads to some uncertainty about the transferability of the identified parameter values for new applications. Furthermore, the model was only implemented in 2D, and the values of C D were taken from experimental observations, which may not always be available. Thus, there is an opportunity to reevaluate the model parametrization, to extend the model into 3D and to incorporate a generic approach to the estimation of C D based on the physical characteristics of the vegetation. Additional work is also required to implement a solute transport model on top of the flow field and turbulence calculations. The following sections review concepts and models relevant to solute transport modeling, before providing an explanation and justification for modifications made to the original King et al. (2012) model parameters.

Solute Transport Modeling in CFD
To model the movement of a solute in a CFD model, an additional scalar transport term and equation is included: where θ is solute concentration and D ij is the diffusion tensor. Equation (12) is analogous to the standard advection-dispersion equation (Fischer et al., 1979). In k-ε models, the diffusion tensor is calculated from the turbulent viscosity as turbulent diffusion: where Sc t is the Schmidt number, the ratio between momentum and mass transport. Sc t is typically assumed to have a value of 0.7 (Tu et al., 2012). D t is applied isotropically due to the turbulence isotropy assumption inherent in k-ε models. White and Nepf (2003) suggested that longitudinal dispersion in vegetation is primarily due to two additive processes: dispersion due to vortex trapping and stem-scale secondary wake dispersion. In vortex trapping, particles are entrained temporarily in eddies formed behind stems. Murphy (2006) simplified the White and Nepf (2003) vortex trapping expression to

Mixing in Vegetation
In secondary wake dispersion, particles travel between areas of lower velocity behind stems and higher velocity between stems, resulting in differential longitudinal advection. Secondary wake dispersion therefore comprises dispersion due to stems (D xxSW ) and dispersion due to gaps (D xxG ). For low-density vegetation (ad < 0.1), White and Nepf (2003) provided a comprehensive approach for determining D xxSW : where σ *2 w = (Γ/16)(C D 3 adR et ) 0.5 , s * = (s + d)/d, s is the distance between stems, R et = (ρUd/(μ t + μ)), and Sc tEFF is the effective turbulent Schmidt number. Γ is an incomplete gamma function, which may be calculated from where c = 2x 0 * C D ad, x 0 * = 0.6, and p = 0.5 (Bhattacharjee, 1970). White and Nepf (2003) also provide a term for the longitudinal dispersion caused by the increase in velocity due to the gap between stems: Tanino and Nepf (2008) described transverse dispersion within vegetation as the sum of turbulent diffusion, and mechanical dispersion, the transverse spread caused by flow path tortuosity: where P snc>r* is the probability the nearest stem is further than r*, where r* is the minimum distance between stems necessary to contribute to lateral mixing, s nc is the center-to-center distance between stems, P snc<5d is 10.1029/2018WR023628

Water Resources Research
the probability the nearest stem is located within 5d, and k ⊥ is the permeability in the direction of flow. Tanino and Nepf (2008) provided expressions for all the parameters in terms of ϕ, assuming a random array of cylinders. This model describes a smooth transition from turbulence-dominated dispersion at lower solid volume fractions to mechanical dispersion-dominated dispersion at higher solid volume fractions.

Objectives
The objectives of the paper are as follows: 1. to implement the turbulence model outlined by King et al. (2012) in 3D; 2. to review the model constants and implement a method for estimating C D from the physical characteristics of the vegetation; 3. to propose a modeling framework suitable for modeling solute transport within vegetation that combines the King et al. (2012) model with the approaches outlined in section 2; and 4. to demonstrate the new model's fitness for purpose using a range of laboratory data sets.

Modeling Framework
This paper proposes a modeling framework that integrates stem-scale mixing processes with the King et al. (2012) k-ε turbulence model to predict solute transport in vegetated flows. The original King et al. (2012) model was developed within a 2D framework, considering a vertical plane, assuming flow in the streamwise direction aligned with the x-axis (orthogonal flow) and assuming that the vegetation could be represented by an array of vertical cylinders of uniform diameter. While it is acknowledged that real vegetation is typically heterogeneous and flexible, arrays of cylinders are regularly used within laboratory studies to represent vegetation (e.g., Nepf, 1999;Serra et al., 2004;Wu & He, 20092009).
In extending the model to 3D, but assuming vertical cylinders, it was assumed that vertical flow does not encounter vegetation, and therefore, the vertical velocity component may be neglected in equations (3), (7), and (8). Stem-scale processes were assumed to occur only within the horizontal, xy-plane. To take into account the nonorthogonality that is possible within a CFD simulation, U was redefined on a cell-by-cell basis for use in equation (7), and so on, as the horizontal, xy-plane streamwise velocity: where u x is longitudinal velocity and u y is transverse velocity. U replaces u i in equation (3).

Drag Coefficient
Assuming vertical vegetation and ignoring skin friction, F z = 0. For drag force in the other directions (equation (3)), C D must be known or estimated. Tanino and Nepf (2008) suggested using an Ergun (1952)-derived relationship to estimate C D in arrays of emergent cylinders representing vegetation: where α 0 is a viscous coefficient, α 1 is an inertial (or pressure) coefficient, and Re d = ρUdμ −1 .  proposed the following parameter values: These values have been found to provide reasonable predictions of C D when ϕ ≤ 0.4 and d ≤ 0.025 m. Equation 21 has been adopted here for emergent arrays. Inertial (or pressure) drag, required in equations (7) and (8), can be obtained from C Dp = 2α 1 . Tanino and Nepf (2008) used the same approach for calculating C Dp , although a different function was used to estimate α 1 .
An alternative approach is required for submerged vegetation. Comparison of results from King et al. (2012) and Tinoco and Cowen (2013), who used the same vegetation in submerged and emergent configurations, respectively, alongside Tang et al. (2014) and Ghisalberti and Nepf (2004), suggests that drag coefficients 10.1029/2018WR023628

Water Resources Research
are reduced for submerged vegetation. Ghisalberti and Nepf (2004) proposed a function relating submerged C Ds to emergent C D , which may be simplified to C Ds = 0.28C D . This ratio has been adopted here to estimate C D for submerged arrays. Tanino and Nepf (2008) showed that turbulent processes within vegetation are affected by vegetation density. Where there is sufficient spacing between the stems, turbulent processes are governed by stem diameter as the relevant length scale. However, once a threshold density of vegetation is reached, stem spacing becomes the relevant length scale. The P w term (equation (7)) assumes that stem diameter is always the dominant length scale when generating turbulence, that is, d < s. While this was correct for the King et al. (2012) test cases, it will not always be so. To form a generic model, equation (7) must be modified to take the length-scale switch into account. King et al. (2012) present the Tanino and Nepf (2008) link between wake-scale turbulence intensity and stem properties as

Length-Scale Switching
where γ is a turbulence-intensity scaling coefficient and by definition γ = (β p /C εD ) 2/3 . Tanino and Nepf (2008), however, took into account the shift in turbulence length scales when s < d, giving where γ 1 = 1.21 and γ 2 = 0.77, with d/s = 0.56 corresponding to the point at which the two curves intersect. d/s = 0.56 is approximately equivalent to ϕ = 0.029 for a random array of cylinders. The first line of this equation is the same as equation (22). Given the relationship between equation (22) and equations (10) and (7), a new equation for P w may be written: where γ ds is the intersection of the two parts of the equation. Equation (24) introduces two new model constants, γ ds and β p2 . Similarly, equation (10) can be reformulated as The modifications to the King et al. (2012) model represented by equations (24) and (25) allow the appropriate turbulence length scale to be applied based on local bulk-scale vegetation characteristics.

Solute Transport Modeling in the Proposed Modeling Framework
As previously stated, when Equation (12) is applied within k-ε models, the diffusion term is the result of isotropic shear-scale turbulence. Applying equation (13) to the King et al. (2012) model, equation (12) therefore includes the effects of both shear-scale and wake-scale turbulence. However, this still does not include the effects of other physical processes (mechanical dispersion, secondary wake dispersion, and vortex trapping) that occur within vegetation when the vegetation is not explicitly modeled at stem scale. The following subsections explain how the expressions given in section 2.4 for estimating dispersion in vegetation are integrated within the proposed modeling framework, resulting in an anisotropic diffusion tensor.

Nonorthogonality
The models for predicting dispersion within vegetation presented in equations (14)- (18) are all nondimensionalized by mean streamwise velocity. Mechanical dispersion is transverse to the direction of flow and of a different order of magnitude to secondary wake dispersion, which is parallel to the direction of flow, and so an anisotropic diffusion tensor must be used. If streamwise flow does not align with the x-direction, that is, flow is nonorthogonal and U ≠ u x , then these equations cannot be applied directly as D ii ≠ D xx and the dispersion coefficients predicted must be rotated to match the orientation of flow.
Following  and equation (19), ψ = tan −1 (u y u x −1 ) describes the angle of the flow within the xy-plane. From this, the diffusion tensor may be written as to take into account nonorthogonality, where and D xx , D yy , and D zz are the x, y, and z components of diffusion with respect to the direction of flow as presented earlier. Kalinowska and Rowinski (2012) and Arega and Sanders (2004) use a similar mechanism to handle anisotropic diffusion within river and coastal flow, respectively. D zz needs no rotation, as it is perpendicular to the xy-plane.

Turbulent Diffusion
To describe turbulent mass diffusion within vegetation, we propose to split equation (11) into separate shearand wake-scale momentum transport equations, consistent with King et al. (2012), giving where μ ts is shear-scale turbulent momentum diffusion and μ tw is wake-scale turbulent momentum diffusion, and μ t = μ ts + μ tw . Equation (13) becomes where D ts is shear-scale diffusion, D tw is wake-scale diffusion, Sc ts is the turbulent Schmidt number where shear-scale turbulence dominates, and Sc tw is the turbulent Schmidt number where wake-scale turbulence dominates. The effective Schmidt number, needed for equation (15), therefore becomes Sc tEFF = (μ t + μ)/(ρ (D ts + D tw )). Choice of Schmidt number will be discussed in section 3.4.
There is currently no explicit model for vertical dispersion within vegetation. Nepf, Sullivan, et al. (1997) and Huang et al. (2008) report vertical dispersion to be less than or equal to transverse dispersion. Lightbody and Nepf (2006) suggest that the ratio of vertical to transverse turbulent diffusion D zz /D yy typically takes values of the order of 0.1-0.2. These observations typically relate to low solid volume fractions, for which overall transverse dispersion is dominated by turbulence-generated dispersion and mechanical dispersion contributes only a minor part. For vertically aligned stems it is reasonable to assume that mechanical dispersion may be ignored. Based on the observed ratios of D zz /D yy of the order of 0.2, it is therefore reasonable to assume that wake-scale turbulence applied vertically may be modeled as 0.2D tw . The effects of shear-scale turbulence are applied isotropically, such that D zz is described by equation (30c).

Complete Mixing Model
Assuming diffusion processes are additive, then the complete dispersion equations used for the anisotropic diffusion within vegetation tensor are as follows: In nonvegetated regions, the diffusion tensor loses the stem-scale processes and becomes

Model Constants
King et al. (2012) suggested values for the five original model constants, β p , β d , C ε5 , C εD , and C λ based on a sensitivity analysis and calibration against laboratory data. King et al. (2012) found their model to be very sensitive to values of C εD and β p and less sensitive to β d and C λ . Model performance was worse when a nonzero value for C ε5 was used, that is, when equation (9) was included. Equation (24) introduces two new model constants γ ds and β p2 . Suitable values for each of these model constants, together with the turbulent Schmidt number, are proposed in this section. The original and modified model constants are presented in Table 1.
The turbulence-intensity scaling coefficient, γ, quantifies the efficiency with which energy dissipated due to stem-induced drag is converted to wake-scale TKE. By definition, γ = (β p /C εD ) 2/3 . The original King et al. (2012) values for β p and C εD imply γ = 0.8. In comparison, Tanino and Nepf (2008) suggested values of 1.21 and 0.77 for γ 1 (d/s < 0.56) and γ 2 (d/s ≥ 0.56), respectively. Both sets of values are based on empirical fits to data. For a generalized model we believe a value of γ = γ 1 = γ 2 = 1.0 provides a reasonable fit to the available data. These values assume perfect efficiency of energy conversion. From this, it follows that γ ds must also take a value of 1.0, as it defines the intersection of the two components of equations (23) and (24). Assuming γ = 1.0 implies that C εD must equal β p . To remain consistent with γ = 1.0 and β d = 1.0, we also set β p = 1.0, assuming perfect conversion of velocity into TKE, and as a result C εD = 1.0 as well.
The constant C λ describes the relationship between wake-scale TKE, wake-scale turbulent dissipation, and turbulent viscosity (momentum diffusion). King et al. (2012) calibrated C λ to a value of 0.01, suggesting that there is little momentum transport due to wake-scale TKE. However, their calibration was not sensitive to this constant. Considering that shear-scale TKE is often significantly greater than wake-scale TKE in shear layers, and their experimental data were dominated by shear layer effects, this insensitivity is expected. Similarly, in uniform vegetation no sensitivity to C λ would be expected, as mean velocities are approximately uniform, and hence, μ t would be constant. A direct comparison of D tw for uniform vegetation with the Tanino and Nepf (2008) turbulent diffusion model suggests that a higher value of C λ is appropriate.
Within the proposed modeling framework, there is no physical process that would suggest wake-scale TKE has a smaller impact on momentum transfer than shear-scale TKE. As such, a value of C λ = C μ = 0.09, as suggested by Launder and Spalding (1974), is proposed.
The last model constants to be discussed are the turbulent Schmidt numbers. The turbulent Schmidt number is decomposed into shear and wake components. A value of Sc ts = 1.0 has been assumed in the unvegetated flow. Preliminary calibration of C λ suggested a minimum value of Sc tw ≈ 0.1 is required. However, the review Table 1 Original and Modified Model Constants

Water Resources Research
of turbulent Schmidt numbers provided by Gualtieri et al. (2017), for a range of environmental flow scenarios, suggests that values less than 0.2 are rarely encountered, and so a value of Sc tw = 0.2 is adopted here. Within vegetation, k s is dampened out, turbulence is dominated by k w , and so this lower-limit Schmidt number is applied when calculating turbuent diffusion. In the unvegetated flow k w is not generated, k s dominates, and so the shear-scale Schmidt number dominates. Near the interface between vegetation and water, k s and k w may be of similar scales, giving an effective Schmidt number between 0.2 and 1.0. This is consistent with the observed value of Sc t = 0.47 reported by Ghisalberti and Nepf (2004) for a shear layer interface.
The modified model constants are therefore β p = β p2 = β d = C εD = γ 1 = γ 2 = γ ds = 1.0, C ε5 = 0.0, C λ = 0.09, Sc tw = 0.2, and Sc ts = 1.0. In the absence of a complete validation data set covering the range of scenarios described by the model, we believe that assuming ideal values for these constants is reasonable. We do, however, test these constants with the modified model against a comprehensive range of experimental data.

Model Implementation
The model has been implemented as a set of user-defined functions in ANSYS Fluent 19 (ANSYS Inc., 2018).
Code for these is available as supporting information to this paper. These functions have been combined with the standard k-ε model in Fluent.
A 5-mm cell size was used with a hexahedral (3D) or quadrilateral (2D) mesh. Vegetated portions of each configuration were designated as porous zones, with C D defined based either on equations (20) and (21) or on laboratory-derived values. Relevant vegetation characteristics (d, s, ϕ, and a) were input within the new user-defined functions.
The free water surface was modeled using a fixed-lid approximation and set as a zero-shear wall boundary. All CFD simulations represent laboratory experiments in a flume and so were constructed as streamwise periodic geometries to allow flow fields to fully develop. All other boundaries were modeled as smooth walls. The enhanced wall treatment function within the CFD package was used (ANSYS Inc., 2018). Depending on cell distance and velocity, this function changes the wall boundary condition from law of the wall to viscous sublayer modeling.
Except where stated otherwise, the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) pressure-velocity coupling scheme was used. The PREssure STaggering Options (PRESTO)! pressure discretization was used with the Quadratic Upstream Interpolation for Convective Kinematics (QUICK) momentum discretization. The turbulence and scalar transport equations used second-order upwind discretization. A relaxation factor of 0.8 was used for the k w transport equation, the same as k s .
All flow fields were first solved as steady state simulations. The flow field was then frozen, and the scalar transport equation solved. For the mixing simulations either continuous or pulse injections were used. To account for the periodicity of the model, an additional dye-removal user-defined function was applied to remove the dye in the 0.1-m length from the upstream and downstream boundaries. Line profiles of concentration were extracted from the CFD simulations at the measurement locations and normalized according to the peak concentration.

Mesh Independence
In a partially vegetated case, strong gradients in the profiles of velocity, turbulence, and dispersion components are expected. It is therefore necessary to confirm mesh independence prior to the presentation of model predictions. Choice of mesh size was verified using a model of the Ghisalberti and Nepf (2005) experiments (see section 4 for detailed analysis of this case), as it has the greatest computational requirements and largest domain of all the models considered here. The domain was meshed using cell sizes of 0.020, 0.010, and 0.005 m with 101,232, 791,198, and 6,429,600 hexahedral cells, respectively. The grid convergence index (GCI; Roache, 1994) was used to evaluate mesh size using a 1.25 factor of safety. Peak velocity, Reynolds stresses, peak shear-scale TKE, and midvegetation-depth wake-scale TKE were evaluated using the GCI.
To evaluate concentration profiles, R t 2 (Sonnenwald et al., 2014;Young et al., 1980) was used to compare the lower-resolution profiles to the highest-resolution profile.

Water Resources Research
largely unaffected (R t 2 > 0.999). Five-mm (or smaller) meshes have therefore been used as they provide sufficient resolution and accuracy for the scale of domains involved in this study.

Illustrative Scenarios
Prior to the specific validation cases considered in section 4, two scenarios are presented here to demonstrate the relative magnitudes of the mixing components that combine to form D xx , D yy , and D zz according to equation (30). The first scenario considers the simplest condition, which is a fully vegetated system, for which the velocity and turbulence properties will be cross-sectionally uniform and the dispersion coefficients will also be spatially uniform. The second scenario introduces a shear layer associated with partial-width vegetation, which in turn leads to spatially varying dispersion coefficients.

Illustrative Scenario 1-Fully Vegetated Case
The main purpose of this scenario is to demonstrate that the current modeling framework (equation (30)) correctly reproduces the original models, and model components, presented by White and Nepf (2003) and Tanino and Nepf (2008) for longitudinal and transverse dispersion, respectively. The original models demonstrate strong dependency on the vegetation solid volume, ϕ. Therefore, simulations have been undertaken to assess how D xx and D yy vary as a function of ϕ. It should be noted that as there is no shear-scale component to the mixing in a fully vegetated case, the D ts component will always be 0.
The models were run as 2D (plan) simulations of a 0.5-m-wide, 0.1-m-long channel, with the enhanced wall treatment and a periodic streamwise boundary. A 1-mm quadrilateral mesh was used. The models were run with d = 0.00635 m, U = 0.030 m/s, Re d ≈ 190, and ϕ varied at 0.01 increments. The mean stem spacing, s, for the vegetation was estimated from ϕ according to Tanino and Nepf (2008), while the frontal facing area, a, was determined from ϕ = adπ/4. These assumptions regarding the estimation of s and a also apply throughout the remainder of the paper. It should be noted that estimates of D xx /Ud are dependent on Reynolds number due to the additional dependence of equation (14) on d and equation (15) on both s and d. Figure 1 presents model predictions of dispersion coefficient as a function of ϕ for a fully vegetated system. Values for the dispersion parameters were taken from the model midchannel; however, the uniform velocity distribution ensures that the dispersion parameters are spatially uniform over all internal mesh elements. Overall differences in D xx between the present model and the original formulations are small. While White and Nepf (2003) suggested that the D xxG term could generally be neglected, it may be seen that its contribution becomes significant for higher values of ϕ, and therefore, it is included in the present model's determination of D xx . The model appears to underestimate D xxSW due to White and Nepf (2003) utilizing fixed values of Re t , C D , and Sc t , as well as their use of an assumption for stem spacing based on a square array, as opposed to a random array of cylinders. The CFD model prediction of D xxVT is exactly the same as in White and Nepf (2003), due to both utilizing equation (14).
Figure 1b similarly confirms that the present model is capable of reproducing the Tanino and Nepf (2008) trends in D yy /Ud as a function of ϕ. As a result of changing γ ds , the peak value of D tw occurs at a higher value of ϕ in the present model, and consequently, it may predict less dispersion at lower stem densities. Regardless, these results confirm that the choice of values for C λ and Sc ts does not significantly affect the magnitude of transverse dispersion.
No plot of D zz has been included here, as equation (30c) implies that its value would simply mimic the pattern of D tw already shown in Figures 1a and 1b. Furthermore, there is no previously reported data on how D zz is expected to vary with ϕ against which to compare the predictions. We acknowledge that this component of the dispersion is subject to greater uncertainties than is D xx or D yy .

Water Resources Research
Comparison of the King et al. (2012) profiles with those generated by the present model confirms that the changes made to the model coefficients (Table 1) do not significantly impact on the estimation of velocity or TKE. The most notable differences are in the calculation of wake-scale TKE, which appears to be dissipated less rapidly with the present model coefficients compared with those proposed by King et al. (2012). This reflects an increased rate of production of wake-scale TKE due to the higher value taken for β p (1.0 versus 0.2) and the increased rate of dissipation caused by the higher value of C ɛD (1.0 versus 0.28). Figure 3 compares the dispersion coefficients produced by the present model with what would otherwise be estimated using a conventional k-ε model. The plots on the right quantify the individual components contributing to the present model. For longitudinal dispersion, two significant differences can be seen. Figure 3b shows the enhanced level of dispersion associated with the vegetated portion of the flow. This primarily reflects the inclusion of secondary wake dispersion (D xxSW ), with lesser contributions due to other processes. Figure 3a shows a spike in the present model predictions of longitudinal dispersion within the vegetation, adjacent to the interface. This spike is a result of increased local contributions, mainly from three processes: vortex trapping (D xxVT ), secondary wake (D xxSW ), and wake-scale (D tw ) dispersion, of around 0.45, 0.30, and 0.20 m 2 /s, respectively. The increases in the predicted mixing coefficients are a result of local increases in streamwise velocity (Figure 2a) and shear-scale TKE (Figure 2b).

Model Validation
The present model has been validated against five published experimental data sets (Table 3). Three types of checks have been carried out to demonstrate that the model is correctly implemented and that it is applicable to a variety of flow conditions. The first set of checks is to confirm that the model correctly predicts the flow field, the second set is to verify that the model correctly predicts solute transport in uniform vegetation, and the third set is to evaluate the model predictions of solute transport in more complex shear flow cases associated with real vegetation.   Velocity and turbulence profiles used to verify model implementation Ghisalberti and Nepf (2005) Artificial random Submerged, vertical shear Validate model on velocity, turbulence, and solute concentration profiles Nepf, Mugnier, et al. (1997) Artificial random Emergent, full width Validate mixing model on D xx Nepf, Sullivan, et al. (1997) Artificial

Water Resources Research
were measured using Nortek acoustic Doppler velocimeters at seven transverse coordinates and six depths. For the submerged vegetation, particle image velocimetry measurements were made in the xz-plane at three transverse positions. Nepf (2004, 2005) provide data suitable for validating the model predictions for flows dominated by vertical shear. They performed experiments on submerged random artificial vegetation in a 0.38-mwide flume with stem diameter 6.4 mm, varying frontal facing area, and flow rate. The vegetation was 0.14 m high, with a constant flow depth of 0.467 m, giving a submergence ratio of 0.3. Velocity profiles were recorded by Ghisalberti and Nepf (2004) using three acoustic Doppler velocimeters separated transversely by 0.1 m where a 0.08-m-long section of cylinders were removed. For the corresponding mixing experiments, a fluorescent dye was released across the width of the channel at the vegetation-water interface, and vertical profiles of concentration were recorded at six downstream locations (0.19, 0.54, 0.92, 1.50, 2.50, and 3.80 m from injection) using digital video cameras.
The experiments of Nepf, Mugnier, et al. (1997) have been used to validate longitudinal mixing predictions under full-width emergent vegetation conditions. Nepf, Mugnier, et al. (1997) placed arrays of random cylinders with a stem diameter of 0.00635 m and varying solid volume fraction in a 0.38-m-wide flume with flow depths between 0.10 and 0.15 m. Dye was continuously injected, with dye concentrations being monitored at two downstream locations using Turner fluorometers at middepth. Longitudinal dispersion coefficients were evaluated from the monitored concentration profiles via solute routing (Fischer et al., 1979).
Similarly, the experiments of Nepf, Sullivan, et al. (1997) were used to validate transverse mixing predictions. Nepf, Sullivan, et al. (1997) placed random arrays of cylinders with different uniform stem diameter and solid volume fractions in a flume. Dye from a continuous injection at middepth and midwidth was monitored downstream using laser-induced fluorescence. Idealized Gaussian profiles of a mass injection (Fischer et al., 1979) were fitted to the observed transverse concentration profiles to obtain transverse dispersion coefficients.
For the transverse shear case, data collected in a partially vegetated channel has been used (West et al., 2018). West (2016) measured transverse profiles of tracer within two types of emergent real vegetation: Typha latifolia (winter and summer). The vegetation was placed to a width of 0.5 m in a 1.0-m-wide channel, and a continuous vertical line source injection of Rhodamine dye was made at the vegetation-water interface. Transverse middepth concentration profiles were measured 1 m and 2 m downstream using a laser-induced fluorescence system. Additional, previously unreported measurements using Carex acutiformis, recorded with the same method, are also presented here. Sonnenwald et al. (2017) characterized this vegetation as follows: winter Typha d = 0.01 m, a = 1.6 m −1 , ϕ = 0.013; summer Typha d = 0.019 m, a = 3.2 m −1 , ϕ = 0.047; and Carex d = 0.005 m, a = 18.3 m −1 , ϕ = 0.077. Stem density was measured by counting the number of stems in a sample area of 0.5 m 2 . This was repeated at up to 10 different locations. At each location, stem diameter was measured at the channel middepth for at least 40 randomly chosen stems, using digital Vernier gauge calipers. As Carex leaves are blade shaped, both width and thickness were recorded. Assuming orientation could be in any direction to the flow, the average of width and thickness was used for diameter. Both live and dead stems were included in the vegetation characterization. Solid volume fraction was estimated as ϕ = π4 −1 Nd 2 and frontal facing area as a = Nd, where N is the stem density. Photographs of the real vegetation test setups are presented in Figure 4, alongside the stem diameter distributions for each type of vegetation, as originally reported in Sonnenwald et al. (2017). Five flow rates (3.4, 4.2, 5.2, 6.4, and 7.5 L/s) were used with a constant flow depth of 0.15 m.

CFD Model Configuration
Details of the CFD model configurations for each of the validation cases are presented in Table 4.
The Nepf, Mugnier, et al. (1997), and Nepf, Sullivan, et al. (1997) models were carried out as 2D horizontal plane simulations 1 m wide by 2.2 m long using a midchannel pulse dye injection of 0.15-s duration with a 0.05-s time step. (A smaller timestep was tested and found to have no impact on the results obtained.) These are judged to be acceptable simplifications; as for emergent simulations the velocity field will be uniform except very close to the walls. Two-dimensional advection-dispersion equation routing, as described by

10.1029/2018WR023628
Water Resources Research Sonnenwald et al. (2017), was used to calculate dispersion coefficients from transverse measurements of dye concentration taken at each time step.
The remaining simulations were carried out in 3D with channel widths and depths as specified in the original laboratory experiments. Channel length was set to allow 0.1 m before the injection and an additional 0.5 m after the furthest downstream measurement point. Whereas in all other cases the SIMPLE pressurevelocity scheme was used, the coupled scheme was used for the Ghisalberti and Nepf (2005) and West (2016) models. Goodness of fit was evaluated using R t 2 . For both sets of model constants, the TKE predictions for the case shown on the left (ϕ = 0.01, d = 3.1 mm) deviate from the experimental data in the upper half of the profile. This is likely to reflect the fact that sparse vegetation (low ϕ) will tend to lead to spatial heterogeneity in the experimental flow field. While the two model predictions are very similar, a small deviation may be observed between the original and modified    Table 1, as well as C D estimated using equations (20) and 21. The quality of the present model prediction is variable over this data set. Typically, the velocity profile is predicted reasonably well, with R t 2 > 0.96 for most cases, although a tendency for velocities to be underestimated in the free water flow and overestimated within the submerged vegetation is noted. This is due to the estimation of C D as a function of velocity, compared with that by King et al. (2012), who used a fixed value. The new formulation of C D reduces the velocity gradient at the interface. This is most evident in the fourth case (ϕ = 0.02, d = 6.2 mm, H/h = 1.89), where the R t 2 falls from 0.997 to 0.942. Neither of the CFD models captures the surface velocity reduction evident at higher submergence ratios, due to simplifications in the modeled free surface boundary condition.

King et al. (2012) Vertical Flow Profile Comparisons
The middle column shows that consistent with the reduced velocity gradient, Reynolds stresses around the vegetation interface are sometimes underpredicted by the present model. TKE (right column) at the interface is similarly underestimated. In most case R t 2 > 0.75, but for the second case (ϕ = 0.02, d = 6.2 mm, H/h = 1.34), the goodness of fit falls to 0.596. However, all predictions of Reynolds stresses and TKE are of the correct order of magnitude, demonstrating that the modified parameter set leads to acceptable predictions of flow and turbulence quantities, even in the absence of a case specific observed C D value.
Figures 6b-6d provide scatterplot comparisons between the experimental and CFD-derived vertical profiles presented in Figure 6a. There are two potential sources for the model deviations noted in Figure 6a: the model constants (see Table 1) and the C D value. The King et al. (2012) model prediction utilized the specific known experimental C D , whereas the present model predictions were based on the generic C D equation presented here as equations (20) and (21). Therefore, a third data set, corresponding to the present model, but using the King et al. (2012) measured C D , is also included in Figures 6b-6d. As was evident in the vertical profile plots, it is clear from the scatterplots that the present model tends to deviate more from the experimental data than does the King et al. (2012) model, underestimating the peak magnitudes of streamwise velocity, Reynolds stresses, and TKE. However, this deviation is almost completely eliminated when the known experimental value of C D is used instead of the present model's estimated value based on equations (20) and (21). The root mean square error values for the data are provided in Figures 6b-6d. The bulk velocity, U A , the quotient of the discharge and cross-sectional area, has been used to nondimensionalize the hydrodynamic parameters in the scatterplots. This comparison confirms that the modified model constants presented in Table 1 do not negatively impact on the present model's capability to predict the flow field and turbulence quantities. The modified model constants have a more notable-and beneficial-impact on the modeling of dispersion, as will be demonstrated subsequently.  Ghisalberti and Nepf (2004) with those generated using the present model. The King et al. (2012) model with our estimated C D values produces identical profiles to the present model, so these have not been included. In the region where streamwise velocity measurements are available, the model predicts a similar shape vertical profile to the measurements. However, as highlighted in Figure 7b (red circles), the present model systematically underpredicts the streamwise velocity. The present model was run using the discharge reported by Ghisalberti and Nepf (2004). This reported discharge is on average 8% lower than the discharge obtained by integrating the velocity profile. To illustrate the impact of this discrepancy on the velocity profile, a mass balanced profile is shown in Figure 7a (blue dot-dashed line). The mass balanced profiles are much closer to the experimental velocities, as evidenced by the improved goodness of fit shown in Figure 7a and the Predicted Reynolds stresses are a good match to peak values at the interface and within the vegetation. However, they overestimate the magnitude as the water surface is approached. Figure 7c shows that most Reynolds stress values are overestimated compared with experimental values. The modeled TKE also demonstrate a good fit to the experimental data (Figure 7d), especially in defining the magnitude of the peak values at the interface at the lower discharges; see Figure 7a (right column). The better overall fits in this case compared with the King et al. (2012) experimental data ( Figure 6) may reflect the higher submergence ratio and overall lower velocities.

Dispersion in Full-Width Emergent Vegetation
Figure 8a compares experimental and CFD-derived longitudinal dispersion for the Nepf, Mugnier, et al. (1997) experiments. It may be seen that the modified model parameters proposed in Table 1 lead to good estimates of D xx at the lower solid volume fractions. For ϕ = 0.055, D xx is overestimated, except at the lowest velocity. This is believed to be due to the high estimated value of C D for these experiments. The model is less sensitive to change in velocity compared with the experimental data.
This comparison highlights inconsistencies between the White and Nepf (2003) and Nepf, Mugnier, et al. (1997) data sets. The model presented here was based on White and Nepf (2003), and it reproduces their observed trends of nondimensional longitudinal dispersion increasing with increasing solid volume fraction (see Figure 1). However, the Nepf, Mugnier, et al. (1997) results show nondimensional longitudinal dispersion reducing with increasing values of both solid volume fraction and velocity. The White and Nepf (2003) laboratory data were analyzed using an increase in variance approach (Fischer et al., 1979), which includes the information provided in the very low concentrations found in the tails of the temporal concentration distributions. While similar experimental configurations were used in Nepf, Mugnier, et al. (1997), here a

10.1029/2018WR023628
Water Resources Research constant injection of tracer was employed. The resulting data were analyzed by fitting a predicted concentration distribution to the first part of the recorded cumulative concentration distribution, as proposed by Chatwin (1980). This approach quantifies the bulk channel flow properties but fails to capture the smaller stem-scale processes, resulting in trends that are contrary to those presented by White and Nepf (2003). Despite these differences, the new model shows good predictions in the midranges of both solid volume fraction and velocity. Figure 8b shows a similar comparison for transverse dispersion with the experiments of Nepf, Sullivan, et al. (1997). Again, the present model replicates the observed data trends reasonably. For 0 < ϕ ≤ 0.017, D yy is estimated particularly well, falling just outside the error bars. For the higher values of ϕ, however, D yy is underestimated. It should be noted that these values are also underpredicted by equation (18) (Sonnenwald et al., 2017;Tanino & Nepf, 2008). Figure 9 presents solute concentration profiles for one of the Ghisalberti and Nepf (2005) submerged vegetation experiments, ϕ = 0.040, Q = 14 L/s. The CFD simulations have used the same solute flux as the experiment, and the results are plotted as relative concentration, with a value of 1.0 equivalent to the peak concentration in the laboratory upstream profile. This experiment has been selected as it exhibits the worst goodness of fit for the present model across the complete data set. At a short distance downstream from the injection, the CFD and experimental concentration profiles match well, with R t 2 values greater than 0.95 over the first 1.5 m downstream of the injection. However, as the distance downstream increases, the CFD concentration predictions diverge from the experimental data, with the model over-predicting the concentration of solute within the vegetation. The poor match at the interface and within the vegetation suggests that shear-driven exchange processes are dominant and that mixing coefficients are underestimated throughout. Overall, the concentration profiles exhibit good agreement with the laboratory data ( Figure 9b). The uncertainty in quantifying vertical dispersion in artificial vegetation, simulated in laboratory experiments by vertical cylinders, is described in section 3.3.2. This contributes to the reduced goodness of fit with distance from the injection point, and the results suggest that the vertical dispersion coefficient within the vegetation D zz may be underestimated. Figure 10 shows the comparison of solute concentration profiles from the real vegetation experiments (West, 2016) with model-generated results. This is the only test case that utilized real vegetation, which accounts for the irregularities in the observed concentration profiles. It should be noted that the upstream plot compares the transverse profiles of dye concentration resulting from a vertical line source injection at the vegetation-open flow interface 1 m upstream, while the downstream predictions presented here were based on the observed-rather than the predicted-upstream concentration profiles. Both upstream and downstream predictions were performed over a 1-m length of channel. Given the inherent variability associated with real vegetation, shown in Figure 4, the modeled predictions are promising. In particular, the extent to which the solute spreads transversely into and within the vegetation appears to be modeled well for the three contrasting vegetation types. Alongside the present model predictions, Figure 10 also shows predictions made using the King et al. (2012) flow field model combined with equation (13) to define dispersion coefficients.

Real Vegetation, Transverse Shear Case
Considering the final case presented in Figure (20) and (21)) was used in both cases. The differences between the two models result from a combination of two interacting factors: the first is the adjustments made to several of the model constants (Table 1); the second is the implementation of the more refined dispersion model (section 3.3). Additional simulations (not shown) indicate that of the two factors, the majority of the observed difference is due to the implementation of the new dispersion modeling framework presented in section 3.3. It should be recognized that for all the simulations of real vegetation, spatially constant vegetation parameters have been employed. This therefore does not account for spatial heterogeneity, patchiness, within natural vegetation.

Discussion
The modeling framework presented here was primarily derived to provide practical estimates of solute transport and mixing within relatively low velocity vegetated aquatic environments, such as ponds and wetlands. Its validity has been demonstrated here for ϕ ≤ 0.1, d ≤ 0.025 m, and Re d ≤ 2,000. Application of the modelusing the parameters proposed in Table 1-in systems characterized by higher stem densities or stem diameters, or with higher flow velocities, should be undertaken only with caution. As many natural ponds include clumps/patches of vegetation, for which the solid volume fraction exceeds 0.1, there is clearly a need for further validation/calibration of the model under these conditions.
The modeling framework also assumes a constant stem diameter, whereas real vegetation is typically characterized by a distribution of stem diameters (Sonnenwald et al., 2017). Details of the vegetation characteristics are provided in Figure 4. A sensitivity test was carried out using one of the summer Typha experiments, varying ϕ, a, and d by ±1 standard deviation. The model appears to be most sensitive to a and d. Frontal facing area, a, influences the drag term, heavily affecting the flow field, while stem diameter, d, influences the dispersion coefficients, affecting solute transport. Interestingly, uncertainty in the solid volume fraction, ϕ, has the least impact. It affects mostly the mechanical dispersion term at higher ϕ values, above the range validated in this study. Similarly, its impact on drag force is greatest at higher ϕ values, and hence, the model has lower sensitivity to low values. Further experimental and theoretical analyses to

10.1029/2018WR023628
Water Resources Research quantify observed and expected differences due to stem diameter distribution effects are therefore strongly recommended.
The modeling approach assumes that the stems behave as rigid cylinders, which are not representative of all types of vegetation. This model has not been validated for leafy vegetation. Use of the leaf area index captures the effects of plant canopies on flow hydraulics (Jalonen et al., 2013). Future work should investigate the use of leaf area index to estimate drag and dispersion. Similarly, the present model does not account for easily deformed, flexible, vegetation or floating wetland systems.
The rationale for the updated parameter values proposed in Table 1 was presented in section 3.4. In several cases, parameter values of 1.0 were proposed, potentially simplifying the model's implementation and usability. The velocity and turbulence profiles presented here have clearly shown that these parameter changes do not negatively impact on the model's flow field predictions.
A notable modification was made to C λ , which was increased significantly from 0.01 to 0.09. While the velocity and turbulence predictions made by King et al. (2012) were insensitive to this parameter, the inclusion of mixing within the present model highlighted strong sensitivities to C λ that necessitated its review and refinement. Preliminary studies using the original King et al. (2012) value of 0.01 (not presented here) revealed poor representation of mixing processes, specifically on underprediction of the extent to which the solute penetrated into the vegetation. A value of C λ = C μ = 0.09, as suggested by Launder and Spalding (1974), was proposed here and has been shown to lead to reasonable estimates of solute mixing.
Following King et al. (2012), the present modeling framework is based on a k-ε turbulence model. All k-ε turbulence models have inherent limitations associated with the assumption of isotropic turbulence. In complex flow fields, such as compound channels, k-ε turbulence models generally fail to reproduce secondary circulations and other important 3D effects. The commercial CFD modeling tool employed here (ANSYS Fluent) includes a full Reynolds stress turbulence model (RSM). In principle, it should be possible to implement the modeling framework outlined in the present paper alongside this or other RSM-based turbulence models.

Conclusions
To the authors' knowledge, the modeling framework presented here is the first CFD k-ε-based modeling approach to accurately represent mixing within vegetation. The framework was developed to provide a practical prediction tool for engineering purposes, and its implementation here within a commercial CFD package potentially facilitates application in complex geometries, such as real, partially vegetated, pond systems.