Computational fluid dynamics modelling of residence times in vegetated stormwater ponds Open Access

Stormwater ponds and wetlands are commonly used in sustainable drainage systems, built to treat rainfall run-off before discharging into watercourses (HA, 2006; Shilton, 2005; Woods Ballard et al., 2015). These systems are rarely monitored post-installation, and sufficient treatment capability on their part is often assumed on the basis of a nominal hydraulic residence time that is described by volume divided by discharge (t_n = VQ⁻¹). A more complete description of the hydraulics is provided by a residence time distribution (RTD), which describes the likelihood a water molecule will stay within the pond for a given duration (Levenspiel, 1972). An RTD reflects the mixing that occurs within a device, and it can indicate short-circuiting, preferential flow paths, recirculation and so on. Figure 1 shows example RTDs for different flow regimes.

An RTD may be obtained experimentally through the use of tracers. Within a pond or wetland, velocities are typically low and tracer duration can be several days (e.g. King, 2006), often making experimental field measurements impractical. An alternative approach is to obtain RTDs through the use of computational fluid dynamics (CFD) modelling. CFD modelling is a flexible tool that allows arbitrary geometries and flow conditions to be represented. Within CFD software, the Navier–Stokes equations are solved using a finite-volume approach over a computational mesh (Tu et al., 2012). Scalar transport modelling or particle tracking allows for an RTD to be calculated (Stovin et al., 2008). CFD modelling has previously been used to evaluate RTDs in a variety of devices, including manholes, stormwater tanks and ponds (Adamsson, 2004; Stovin et al., 2013; Tsavdaris et al., 2013).

Stormwater ponds and wetlands contain vegetation, which can significantly affect flow hydrodynamics (Nepf, 1999). The most obvious impacts of vegetation are that it acts as an obstacle within the flow and that the flow velocity is significantly reduced within vegetated zones. Saggiori (2010) and Tsavdaris et al. (2013) proposed that vegetation within a pond could be modelled within CFD by representing it as a porous zone, where drag is calculated based on the stem diameter (d) and the solid volume fraction (ϕ) using the Ergun equation (Ansys, 2014). Saggiori (2010) used the porous zone approach to successfully reproduce velocity profiles around submerged vegetation in a channel, as well as conducting an initial study into producing RTDs within vegetated flows using a Lagrangian particle tracking approach. Tsavdaris et al. (2013) confirmed that the porous zone model reflected large-scale hydrodynamics within a pond around a patch of vegetation.

However, in addition to its control over the large-scale flow field, vegetation also generates stem-scale turbulence that leads to enhanced mixing within the vegetation. Transverse mixing processes at the stem scale include turbulent diffusion and mechanical dispersion (Nepf, 1999), while longitudinal processes include vortex trapping and secondary wake dispersion (differential shear) around stems (White and Nepf, 2003). Stem-scale mixing processes are not represented within the porous zone model. Sonnenwald et al. (2016a) showed that the porous zone approach reproduced solute traces within channels containing submerged vegetation (where the dominant mixing effects are due to the large-scale geometric control imposed by the vegetation/clear water boundary), but not within emergent vegetation (where stem-scale mixing becomes dominant). The k–ε turbulence modelling approach only reproduces larger scale dispersion effects due to velocity shear. As the velocity profile within emergent vegetation is near uniform in the direction of flow (Nepf et al., 1997), there is no velocity shear and therefore no mixing is modelled.

King et al. (2012) proposed a k–ε turbulence model for vegetation, which will be referred to hereafter as the KTC model. The KTC model is based on the standard k–ε model, incorporating a second stem-scale turbulent kinetic energy term that uses vegetation characteristics to estimate turbulence occurring within the vegetation. As it introduces turbulent kinetic energy within the vegetation, it should therefore also account for mixing within the vegetation. Sonnenwald et al. (2016b) compared the KTC model with a generalised form of the porous zone model referred to as the momentum sink approach, where the Ergun-derived parameters of the porous zone approach are replaced with a bulk drag coefficient (C_D), frontal facing area (a) and solid volume fraction modification (King et al., 2012; Nepf, 1999). Sonnenwald et al. (2016b) showed that the KTC model resulted in enhanced levels of mixing within the vegetation compared with the basic momentum sink (or porous zone) representation. However, the KTC model only represents the turbulent dispersion process described by Tanino and Nepf (2008), applied isotropically. The KTC model does not represent transverse mechanical dispersion, longitudinal vortex trapping or longitudinal secondary wake dispersion, all of which are relevant to mixing within vegetation (Sonnenwald et al., 2017).

This paper describes an alternative, empirical, approach to incorporating the smaller scale vegetation-driven dispersion processes within a CFD model. Sonnenwald et al. (2017) recorded experimental values of dispersion within real, emergent, Typha latifolia vegetation. In this paper, the momentum sink approach is combined with these experimental values to produce RTDs using CFD simulations for four typical stormwater pond and wetland designs, both with and without the inclusion of vegetation.

The commercial CFD package Fluent (Ansys, 2014) was used to model four vegetated stormwater ponds in three dimensions: a regular basin (RB), a curvilinear basin (CB), a pond with island (PI) and a wetland (W). These ponds, shown in Figure 2, were based on design guidance provided by the Highways Agency (now Highways England) (HA, 2006). All ponds had a surface area of approximately 500 m² with a 2·0 m wide by 0·1 m deep inlet and 1·5 m wide by 0·15 m deep outlet. Side slopes were 1:3. The pond depth was 1·5 m, while wetland depth alternated between 0·45 and 0·60 m every 1 m along its length. Wetland depth transitions also had a 1:3 slope. The models were assumed to be in a ‘storm’ scenario, with a flow rate of 30 l/s and an inlet velocity of 0·15 m/s.

Two different vegetation growth stages were used, a winter and summer T. latifolia, described in Table 1 (Sonnenwald et al., 2017) and shown in Figure 3. The summer vegetation, while having a similar stem density to the winter vegetation, had larger stem diameters and a solid volume fraction three to four times greater. Vegetation within the ponds was positioned on the sloping sides. As inlet configuration can have a large impact on the flow field, a reduction in plant density near the inlet (e.g. due to scour) was not simulated; instead, the vegetation was of uniform density. The entire wetland model was vegetated and, additionally, had two baffles (1·0 × 0·5 m in plan) near the inlet. It is worth noting that the regular basin pond and the curvilinear pond were both characterised by a central clear flow channel that ran almost from the inlet to the outlet; the pond with the island and the wetland were characterised by cross-sections where the vegetation extended across the full width. One set of vegetated ponds (with winter T. latifolia) was modelled with stem-scale mixing omitted. This provided a comparison between the previous porous zone models (without stem-scale mixing) and the current approach. The models were also run without vegetation, giving a total of 16 models.

Fluent (Ansys, 2014) simulates fluid flows using an iterative finite-volume approach to solving the Navier–Stokes equations. These are the continuity equation

and the momentum equation

where ρ is density, t is time, u_i is instantaneous velocity, p is pressure, $Tij$ is the stress tensor, F_i represents any external force acting on the control volume and i (or j) represents the three dimensions x, y and z in tensor notation (Tu et al., 2012). In the case of the momentum sink approach, the force of the vegetation acting on the flow is described by

where the over-bar indicates local mean velocity (Sonnenwald et al., 2016b). The Fluent porous zone model (not shown here) is analagous to Equation 3 when the porous zone viscous drag coefficient α is zero (which disables the viscous term) and the porous zone inertial drag coefficient C₂ is defined by

Drag force is applied only within the vegetation and is assumed equal in the x- and y-directions and zero vertically.

Tanino and Nepf (2008) and Tinoco and Cowen (2013) experimentally evaluated the drag force due to cylinders representing vegetation and presented empirical results relating vegetation characterisation to drag. For each vegetation tested they obtained best-fit values for α₀ and α₁ in the relationship

where $Red=Udν−1$ is the stem Reynolds number, U is the mean streamwise velocity and ν is the kinematic viscosity of water.

Sonnenwald et al. (2017) combined the results obtained by Tanino and Nepf (2008) and Tinoco and Cowen (2013) to provide estimates of α₀ and α₁ based on vegetation characteristics

To estimate the drag coefficient, α₀ and α₁ from Equation 6 were substituted into Equation 5. C_D was limited to a maximum value of 10 for numerical stability.

A virtual dye injection was carried out to obtain an RTD from each of the CFD flow field simulations. The virtual dye was modelled using the scalar transport equation

where θ is concentration, μ_t is the turbulent viscosity (the rate of momentum transport due to turbulent stresses) and Sc_t is the turbulent Schmidt number (the ratio between momentum and mass transport). This is directly analogous to the advection–diffusion equation (Rutherford, 1994), where the diffusion coefficient is the ratio of the turbulent viscosity to the Schmidt number (⁠$D=μtSct−1$⁠).

The new approach presented here modifies the scalar transport equation to include an anisotropic diffusion tensor describing dispersion within the vegetation, giving

where $Dij$ is the new stem-scale vegetation dispersion term. $Dij$ is zero in clear water.

Dispersion due to vegetation is dependent on the mean streamwise velocity (Sonnenwald et al., 2017). As such, linear relationships between U and longitudinal dispersion (D_xx), transverse dispersion (D_yy) and vertical dispersion (D_zz) need to be incorporated into $Dij$⁠. A few simplifying assumptions are necessary: first that the vegetation is vertical and second, as a result, it is primarily flow acting perpendicular to the vegetation (i.e. in the horizontal plane) that causes dispersion due to vegetation. Following on from this, $Dij$ is given by

If the mean streamwise velocity is in the x-direction, k₁ is D_xx, k₂ is D_yy and k_ψ is 0. However, a pond may be any geometry and therefore contain significant non-orthogonal flow. As streamwise velocity may not be aligned with the coordinate system, it is necessary to calculate the shear components of the longitudinal and transverse dispersion. Local streamwise flow velocity may be calculated as $uψ=(u¯2+v¯2)1/2$⁠, where v is the instantaneous transverse velocity and $ψ=tan−1⁡(v¯u¯−1)$ describes the angle of the flow with respect to the coordinate system. Therefore, after applying a rotation matrix, the horizontal plane diffusion components become

where D_ii = B_iiu_ψ and B_ii is the slope of the linear relationship between velocity and dispersion taken from Table 2. Sonnenwald et al. (2017) did not make direct measurements of D_zz, so an estimate of B_zz was made as B_yy/3 assuming D_zz < D_yy < D_xx, based on the work of Nepf et al. (1997).

Model geometries were meshed with 0·02 m hexahedral cells for the inlet and outlet channels. The main body of each model was meshed using tetrahedral cells and the advanced sizing function on proximity and curvature, with a minimum cell size of 0·05 m and a maximum cell size of 0·3 m. This placed more cells where the geometry was more complex (Ansys, 2014). These cell sizes were chosen to ensure that all the main geometric elements (e.g. inlet width or pond depth) were described by a minimum of five cells. The models each have of the order of one million cells. One of the simulations was run at an increased mesh resolution to test mesh independence. All of the major flow features were reproduced, although local velocity magnitudes changed on the order of 5%. Although full mesh independence would be desirable, the higher resolution was practically infeasible due to the availability of computational resources.

Second-order discretisation was used along with the k–ε realisable turbulence model (Tu et al., 2012). This realisable model was chosen as it reliably reproduces jet-driven recirculation (Stovin et al., 2013). The enhanced wall treatment was used to account for both the higher velocities near the inlet and low velocities within the body of the pond. A fixed-lid approximation was used in which the free water surface was modelled as a symmetry boundary. After solving a steady-state flow field, a pulse injection of dye was traced through the stormwater pond using a time-dependent simulation. The pulse injection was made by setting the dye fraction in the inlet channel to 1, leaving it at 0 for the rest of the pond. The RTD was then determined by recording the average dye fraction across the outlet face at 5 s intervals. As a pulse injection was used, the temporal concentration profile at the outlet directly represents RTD after being normalised to have an area of 1. The simulation was run with a time step of 0·1 s for a simulation duration of 1·5 d. The turbulent Schmidt number was assumed to be 1.

An RTD is unique to each system's geometry and flow rate. Danckwerts (1953) showed that by normalising time by the nominal residence time (i.e. the normalised time is t/t_n) RTDs representing different systems could be compared independently of flow rate and volume.

While the RTD or cumulative RTD (CRTD) provides a complete description of a system's mixing characteristics, a number of metrics – derived from the CRTD – are commonly used to quantify system hydraulic performance. Teixeira and do Nascimento Siqueira (2008) and Farjood et al. (2015) reviewed several of these metrics. In general, there are two types of hydraulic performance index, those describing short-circuiting and those describing mixing. An ‘ideal’ treatment system might be considered to have no short-circuiting and little mixing (close to plug flow) such that all pollutants experience a consistent residence time and are therefore assumed to experience a consistent level of treatment (Persson et al., 1999).

The choice of appropriate hydraulic performance indices is a complex subject and not the focus of this paper, so only fundamental parameters of the CRTD are considered here. λ_x is used to denote times normalised by t_n (i.e. λ_x = t_x/t_n). λ₁₀ (t₁₀/t_n) is the normalised travel time for the first 10% of the dye (i.e. the normalised time taken for the first 10% of the dye to pass through the system). λ₁₀, λ₅₀, and λ₉₀ may be used to describe the CRTD. The closer each of these values is to 1, the closer the flow field is to plug flow. Low values of λ₁₀ and λ₅₀ combined with high values of λ₉₀ typically indicate the presence of short-circuiting and dead-zones.

Figure 4 shows velocity contours and vectors for the ponds with no vegetation. Relatively strong jets are formed in the regular basin pond, the pond with island and the wetland. The CFD simulations for these devices do not quite resolve to steady state, with eddy shedding occurring within the recirculation zone (not shown). The inlet jet within the curvilinear basin pond is much shorter, due to it impacting on the side wall. The reflection of the jet from the sides of the pond allows a more stable flow field to develop. The establishment of a single recirculation cell in the empty wetland (effectively just a basin), despite its symmetric geometry, is expected from the asymmetric meshing and inlet jet oscillation (due to high jet Reynolds number). Adamsson et al. (2003) reported a similar asymmetric flow field in a symmetric detention tank model.

Figures 5 and 6 show velocity contours and vectors within the stormwater ponds vegetated by winter and summer T. latifolia, respectively. Vegetation surrounded the edges of the ponds and was uniform throughout the wetland. The flow fields for both seasons of vegetation are very similar, suggesting that geometry plays a larger role than vegetation characteristics in determining the flow field. However, the increased density of the summer T. latifolia is reflected in the visibly increased momentum diffusion around the inlets (larger patches of the same velocity magnitude).

Figure 7 shows RTDs for the ponds without vegetation. All four devices show large peaks in the RTD at t/t_n < 0·2, which is well below the nominal residence time. The highest short-circuiting peak, observed in the RTD of the pond with an island, corresponds to the longest inlet jet visible in the flow field shown in Figure 4(c). For the four devices considered here, the magnitude of the RTD peak varied approximately in proportion to the length of the inlet jet.

Figure 8 shows RTDs for the vegetated devices, while Figure 9 provides a comparison between the CRTDs for non-vegetated cases and the summer T. latifolia. Both types of vegetation produced very similar RTDs. For the three ponds, the peak of the RTD with winter T. latifolia is slightly greater than summer T. latifolia. This is reversed for the wetland. All peaks in the RTDs for vegetated systems occur at t/t_n > 0·3, indicating less short-circuiting than the non-vegetated devices. The peaks are also of lower magnitude. The regular basin and curvilinear basin ponds performed similarly once vegetated. In both cases, the first arrival times and RTD peaks are later than when unvegetated, but earlier than for the vegetated pond with island or the wetland. The wetland RTD is centred around t/t_n = 1 and shows the smallest spread, indicating that it is close to idealised plug flow.

The RTD metrics are presented in Table 3 and Figure 10. These are consistent with the flow field and RTD comparisons presented above. All λ₁₀, λ₅₀, and λ₉₀ values move closer to 1·0 with the inclusion of vegetation, indicating closer to plug flow. The increase in λ₅₀ is minimal for the curvilinear basin. The low λ₅₀ value of 0·06 for the unvegetated pond with island reflects the extreme short-circuiting of this configuration. The small decrease in λ₉₀ for the regular and curvilinear basins indicates the continued presence of dead-zones.

All the RTDs for the non-vegetated stormwater ponds have long tails, which, although not shown in Figure 9, are indicated by the λ₉₀ values in Table 3 and Figure 10. The RTDs reflect the occurrence of dead zones, the centre of the recirculation cells and, in some cases, the corners of the devices. Tsavdaris et al. (2014) showed similar recirculation zones in non-vegetated ponds. German et al. (2005) modelled a non-vegetated stormwater pond with an island and showed a recirculation zone similar to the one visible in Figure 4(c) formed around their island. Interestingly, removing their island did not alter the RTD. The effect of recirculation is also visible in the secondary peaks in the RTDs. Of the non-vegetated devices, the curvilinear pond exhibited the least short-circuiting. Its performance is similar to the idealised cells-in-series model (Levenspiel, 1972) and reflects the daisy-chaining of the three recirculation cells that is established between the bends of the pond.

In general, the presence of vegetation led to more uniform flow, as is evident in Figures 5 and 6. These models did not form recirculation zones and this is reflected in the absence of secondary peaks in their RTDs. Vegetation had the least impact in the case of the curvilinear basin pond, suggesting that the complicated ‘S’ shape of the geometry has a strong influence on the flow field. The conceptual design guidance behind this basin and the regular basin is that there is a 1:3 width to length aspect ratio to ensure sufficient flow path length for treatment, regardless of geometry (Woods Ballard et al., 2015). However, it is evident that the presence of vegetation is critical to achieving this in the regular basin pond. The contrasting RTDs shown for the wetland in Figures 7 and 8 suggest that the presence of vegetation is also critical to the performance of the wetland. Persson et al. (1999) showed a similar increase in hydraulic efficiency when a wetland was uniformly vegetated compared with a base case of sparse vegetation near the edges. The vegetated pond with an island performed similarly to the wetland. In this case, the island caused the vegetation to extend across the entire width of the pond, ensuring that there was no higher velocity flow path in the clear water areas between the inlet and outlet.

It is evident that the effects of geometry (i.e. pond shape and vegetation layout) are fundamental to determining flow fields and the corresponding RTDs. These observations are consistent with several previous studies (e.g. German et al., 2005; Khan et al., 2013). However, while previous researchers did not attempt to model stem-scale mixing processes occurring within the vegetation, these effects were represented in the current work. Differences between winter and summer T. latifolia are readily apparent in the increased momentum diffusion within the near-inlet flow fields of the ponds (Figures 5 and 6). The relative importance of including these stem-scale processes is highlighted in Figure 11, which contrasts RTDs generated from models including stem-scale effects with RTDs generated from models that do not include the addition of stem-scale mixing processes. The latter were modelled using the momentum sink (or porous zone) approach without the additional anisotropic dispersion component in Equation 8. It may be seen that the RTDs corresponding to the models that do not include stem-scale processes are very similar to those that do. For the regular and curvilinear basins they are nearly identical while, for the pond with an island and the wetland, the exclusion of stem-scale mixing processes results in a greater peak concentration (10% and 20%, respectively) and a slightly later first arrival time. The inclusion of stem-scale mixing processes can have a significant impact on the overall performance of the device, depending on vegetation layout. While the similarities between Figures 5 and 6 and Figures 7 and 8 suggest that the specific characteristics of the vegetation may be less critical than whether or not it is modelled at all, the data demonstrate, for the first time, that the effects due to mixing within the vegetation also need to be properly represented, particularly in devices that include full-width vegetation.

The current approach to modelling mixing within the vegetation is empirical. The dispersion parameters applied in this study were for one specific type of vegetation (T. latifolia) and the values used for D_zz were estimated, not measured. The vegetation was assumed to be perfectly vertical, with dispersion primarily occurring in the horizontal plane. There was also a sharp boundary between the clear water and the vegetation; in reality, this would be a rougher interface generating additional turbulence and local flow effects. Velocity shear around this interface appears to be the dominant source of mixing within the current models, so this is a critical aspect to address. Further modification of the KTC model (King et al., 2012) to include the stem-scale processes described empirically here may facilitate a more complete description of the mixing both within the vegetation and at the interface between the vegetation and the free water.

In addition to these theoretical limitations, the assumption that flow at the inlet would pass directly into vegetation is unrealistic as scour will typically prevent vegetation from growing. This assumption strongly influences the flow field, leading to the diffusion of inlet momentum and a reduction in potential short-circuiting flow between the inlet and outlet. In contrast, Tsavdaris et al. (2014) modelled a clear inlet leading towards vegetation and showed increased turbulence intensity due to the vegetation obstructing the flow. Due to this sensitivity, additional investigations of inlet conditions would be of practical benefit.

In addition, the models only considered a high flow rate associated with storm conditions whereas, in reality, conditions vary significantly as a pond fills and empties and not all stormwater ponds have base flows. However, time-varying CFD simulations are both complex to configure (especially, for example, when considering changes in water level) and computationally expensive. Werner and Kadlec (1996) suggest that it may be possible to adapt steady-state RTDs to these more complex time-varying scenarios and this idea warrants further exploration.

All of the vegetated stormwater ponds exhibited satisfactory performance, with λ₅₀ > 0·5. The λ₅₀ values of the pond with the island and the wetland were higher, with the wetland approaching λ₅₀= 1. λ₅₀ values of around 0·5 suggest that flow is only experiencing half the available volume and therefore total volumes should be up to 100% larger to accommodate this when designing against treatment criteria. All of the non-vegetated systems performed noticeably less well than their vegetated counterparts, although the curvilinear pond performed well even without vegetation due to a combination of its long flow path and interconnected small-scale recirculation zones. Persson et al. (1999) showed a non-vegetated rectangular pond with flow directed by baffles into a sinusoidal shape to be among the best performing pond layouts. Wols et al. (2010) and Chang et al. (2016) also reached similar conclusions regarding baffles or sinusoidal flow paths. The use of curvilinear pond outlines and/or baffles may therefore mitigate the risk of a system under-performing due to uncertainties in the quality and seasonality of the vegetation.

A new approach to modelling mixing due to vegetation within CFD models was developed, enabling stem-scale mixing effects occurring within a patch of vegetation to be represented for the first time. Four stormwater ponds were modelled based on guidance issued from Highways England (a regular basin, a curvilinear basin, a pond with an island and a uniformly vegetated wetland) to obtain RTDs for estimating hydraulic efficiency. The stormwater ponds were modelled with two seasons of one type of vegetation (winter and summer T. latifolia) as well as without vegetation. Accounting for the presence of vegetation within the CFD models of the ponds and wetlands significantly affected the simulated flow fields and, as a result, their RTDs. All the non-vegetated ponds showed significant short-circuiting compared with vegetated ponds. The new modelling approach reflected the effects of season, although the two seasons of vegetation gave similar performance. It is argued that while it is critical to include vegetation within a model geometry and to correctly account for mixing processes within the vegetation, accurate characterisation of the vegetation may be less critical. Further refinement of the modelling approach is recommended to address the simplifications currently applied to the representation of the water/vegetation interface.

Based on the new modelling approach, the current design guidance appears to be sufficient for producing functional stormwater treatment ponds and wetlands, but volumes should be between 50% and 100% larger than those suggested by nominal residence times where the primary design criterion is stormwater treatment.

This research was funded by the Engineering and Physical Sciences Research Council (grants EP/K024442/1 and EP/K025589/1) and the Warwick Impact Fund. The authors would like to thank colleagues at Highways England for their insights into pond design.

a

frontal facing area

B

dispersion slope

C₂

inertial drag coefficient

C_D

bulk drag coefficient

D

diffusion coefficient

$D$

diffusion tensor

d

stem diameter

F

external force acting on the flow

k

non-orthogonal diffusion coefficient

N

stem density

p

pressure

Q

discharge

Re_d

stem Reynolds number

Sc_t

turbulent Schmidt number

t

time

t_n

nominal residence time

t_x

time taken for x% of dye to pass outlet

U

mean streamwise velocity

$u¯$

mean instantaneous velocity

u, v, w

instantaneous velocity

V

volume

α

viscous drag coefficient

α₀, α₁

best-fit coefficients

λ_x

t_x in normalised time

θ

concentration of solute

μ_t

turbulent viscosity

ν

kinematic viscosity

ρ

density

$T$

stress tensor

ϕ

solid volume fraction

ψ

streamwise flow angle