An Experimental Investigation of the Hydraulics and Pollutant Dispersion Characteristics of a Model Beaver Dam

Abstract

Beavers have influenced the world’s ecosystem for millions of years. Their dams create ponds and wetlands that provide a large range of hydraulic and ecological benefits to the natural world, including mitigation against flooding and improving water quality. As beavers are now being reintroduced to many parts of the world, it is important to fully understand the impact of their dams on the flow characteristics of the water-courses on which they are built. This paper investigates the relationship between the physical properties of a model beaver dam and its fundamental hydraulics and pollutant dispersion characteristics. The first objective of this paper was to develop a modelling framework to relate discharge to flow-depth for dams with a combination of porous and impermeable sections. The second objective was to utilize a similar framework to predict the down-stream concentration distribution of an up-stream pollution event passing through such systems. The ability to model these parameters for dams with variable lengths of porous and impermeable sections is important as the porosity of beaver dams can vary with depth, depending on which sections are constructed from branches, rocks, or compacted mud. The analysis and modelling developed in this paper show that a single, general relationship can be obtained between discharge and flow-depth regardless of the presence of sections that are both porous or impermeable, provided the relative depths of these sections are known and accounted for. It is also shown that the Nominal Residence Time and the Advection Dispersion Equation can be used to predict pollutant transport in such systems. These two equations have previously been shown to have limitations when applied to some complex systems, so demonstrating they can be applied to a porous dam with combinations of porous and impermeable sections at the relative discharges investigated is noteworthy.

1. Introduction and Literature Review

1.1. Background

Beavers have played a major role in shaping the world’s ecosystem for millions of years. A remarkable and unique attribute of the species is their ability to construct dams from tree trunks, branches, mud, and rocks [1]. These dams slow stream velocities and create large areas of storage in the form of beaver ponds and wetlands [1,2,3,4,5,6].

Beavers primarily build dams to create a safe habitat from predators and to allow for easy and safe access to food [7]. However, their ponds also have an enormous positive secondary effect on the surrounding ecosystem. Wetlands created by beaver ponds provide habitats for a whole range of species [1]. Beaver dams also perform key hydrologic functions such as dissipating stream energy, attenuating peak flows, and increasing groundwater recharge and retention [8,9]. These increase summer low flows and elevate groundwater levels in stream valleys, thus expanding the extent of riparian vegetation [10,11]. For example, in the headwaters of the Colorado River in Rocky Mountains National Park, beaver dams have increased surface water and groundwater in both high-low flow periods, attenuating declines in the water table during drier periods [11,12].

Beaver Dams also have a positive impact on water quality. They can aid in the treatment of pollutants through a combination of natural filtration and increased transient storage and residence times [13,14,15,16,17,18,19,20,21,22]. Quantifying transient storage and solute transport in natural watersheds is of growing interest when addressing the effects of climate change [23]. It is therefore crucial to establish relationships that can link dam porosity, solute transport and storage to flow conditions.

The positive impact of beaver dams has led to the use of human-made ‘leaky barriers’, also referred to as woody debris, woody debris structures, woody debris jam, and leaky dams [24]. Leaky barriers are constructed out of large logs and use the principle of ‘biomimicry’, replicating the effect of natural woody material accumulation in streams, including windfall, snagging by hanging riparian vegetation, and beaver dam construction. However, unlike beaver dams, there is often no ‘on-site’ maintenance for leaky barriers, and as such, they are constructed to allow low-flows to continue unimpeded to prevent impoundment of water and only enable attenuation and increased flow residence during spate conditions.

1.2. Previous Work on Beaver Dams and Leaky Barriers

Previous work on beaver dams has investigated sedimentation processes in single beaver ponds [25,26,27,28,29] or with sequences of dams [13,30,31]. However, the specific effects of a beaver dam depend primarily on the magnitude of flow. Despite the current understanding that little geomorphic change occurs during base flows [32], past research is limited to data collected under such conditions. The closest study to flood discharges was Daniels and Rhoads [33], who studied the three-dimensional flow structures around large woody debris for two flow stages. However, the higher discharge was below bankfull conditions, and so there is still a lack of empirical documentation of how beaver dams affect flow at bankfull or flood stages.

A number of 1D methods have been used to quantify the hydraulic and hydrological effect of leaky barriers. A simple method often employed is to establish a value for the Manning’s roughness value ‘n’ based on the leaky barrier [34,35]. However, these values cannot easily be generalized, as they are often specific to site and flow conditions induced by specific rainfall events [34,36]. Weir equations can also be used as a simple approach to estimate storage, by representing leaky barriers as different shaped weirs [24,37] or as sluice gates [24,38,39]. However, the sparsity of field data available makes it challenging to correctly identify datasets to be used for the estimation of boundary conditions to calibrate and validate numerical models, and there is a growing need to design experimental facilities that can replicate flow through leaky barriers [24]. Leakey et al. [24] developed a 1D Godunov-type scheme to model leaky barriers that compared well with laboratory data.

In a similar fashion to beaver dams, leaky barriers also impact pollution transport and dispersion due to the creation of large storage ponds [40,41,42,43,44]. Leaky barriers generate temporary retention and gradual release of solutes, creating an asymmetric shape in the concentration vs. time curves, leading to the possibility of reactive pollutants being absorbed by sediment [45]. Additionally, multiple studies have confirmed that these structures have a significant effect on transport behavior, especially in mountain streams [22,46,47,48,49].

1.3. Pollution Transport Modelling

A useful starting point for modelling pollution transport in any system, especially where storage is a primary consideration, is to predict the system’s residence time i.e., the average amount of time a solute will ‘reside’ in the system. An estimate of residence time can be made using the ‘Nominal Residence Time’ (NRT), which assumes ideal plug-flow through the system:

$$\mathrm{NRT}=\frac{\mathrm{V}}{\mathrm{Q}}$$

(1)

where NRT is the Nominal Residence Time (s), V is the system’s volume (m³), and Q is the discharge (m³/s).

For systems that approximate ideal plug-flow conditions, the NRT gives a simple and robust estimate of the residence time. However, recent studies have shown that the NRT provides a poor estimate for more complex systems [50,51,52,53,54]. For example, in systems with overgrown vegetation, ‘short-circuiting’ caused by preferential flow paths around the dense vegetation can significantly reduce the working volume of the system and lead to residence times that are far lower than the NRT [55].

Predicting the down-stream concentration distribution of solutes for a pollution event can be achieved using the modelling framework of Taylor [56,57]. Taylor’s classical work showed that the down-stream concentration distribution of a solute, introduced at some point up-stream, can be predicted using an adapted version of Fick’s second law of diffusion. Here, Fick’s molecular diffusion coefficient is substituted with the ‘Longitudinal Dispersion Coefficient’, a coefficient that determines how much longitudinal spread has occurred from the up- to the down-stream location. Taylor’s standard solution to the Advection Dispersion Equation is given below (deemed ‘ADE framework’ in this paper). The ADE framework allows for the prediction of a down-stream concentration profile based on an estimate of the system’s residence time and the Longitudinal Dispersion Coefficient. This provides a 1D, cross-sectional average concentration vs. time profile that accounts for the dispersion effects that have taken place over a defined reach.

$$c\left(x,t\right)=\frac{M}{A\sqrt{4\pi D_{xx}\overline{t}}}\left[\frac{{(x-ut)}^2}{4D_{xx}t}\right]$$

(2)

where $c\left(x,t\right)$ is the concentration profile prediction (ppb) at location $x$ (m) and time$t$(s), $M$ is the mass of the solute (kg), $A$ is the flow’s cross-sectional area (m²), $D_{xx}$ is the Longitudinal Dispersion Coefficient (m²/s), and $u$is the mean velocity (m/s).

Equation (2) assumes an initial ideal ‘slug’ distribution of the solute and predicts a down-stream profile that is perfectly Gaussian. For a non-ideal, known up-stream pollution event, Equation (2) can be modified to predict a down-stream profile based on the known up-stream distribution [58]. Here the down-stream profile is not necessarily Gaussian, but the transfer function between the up- and the down-stream profile remains Gaussian:

$$c\left(x_2,t\right)={{\displaystyle \int }}_{\gamma -\infty }^{\infty }\frac{c\left(x_1,\gamma \right)u}{\sqrt{4\pi D_{xx}\overline{t}}}\exp \left[\frac{u^2{\left(\overline{t}-t+\gamma \right)}^2}{4D_{xx}\overline{t}}\right]\mathrm{d}\gamma$$

(3)

where $c\left(x_2,t\right)$ and $c\left(x_1,\gamma \right)$are the up- and down-stream profiles, respectively (ppb); $u$is the mean velocity (m/s); $t$ is time (s); $\gamma$ is an integration variable (s); $D_{xx}$is the Longitudinal Dispersion Coefficient (m²/s); and $\overline{t}$ is the travel time (s).

This classical model has been shown to work well under many conditions, but similar to the NRT, the ADE framework has also been shown to have limitations for more complex systems and flows, mainly due to the assumption of a Gaussian transfer function [59,60].

1.4. Summary

It is clear from the literature that further work is required to establish a simple and general relationship between the physical structure of beaver dams (and analogous structures such as leaky barriers) and the systems underlying fundamental hydraulic parameters, such as discharge and corresponding up-stream flow-depth. These parameters are key to understanding the positive or negative effect of such structures on flooding. Furthermore, there is also a need to link the physical characteristics of the dam to its effect on pollutant transport and dispersion, to be able to assess how such structures can affect water quality. It is also important to be able to understand these parameters in the context of a dam with both porous and impermeable sections, as is the case for beaver dams where sections of the dam made out of compacted mud and rock can often be almost impermeable [1].

The aim of this paper was to experimentally investigate the underlying hydraulics and pollution transport characteristics of a model porous dam designed to simulate the flow behaviour of a beaver dam or a leaky barrier. The data collected was used to establish a semi-empirical model that, at a given flow-depth, can predict the system’s discharge, the dam’s pollutant residence time, and the pollutants’ down-stream concentration distribution for a given up-stream pollution event. A relationship between discharge and flow-depth is established that is generalized for dams with any combination of porous or impermeable sections. This relationship is extremely important to understand how such systems would respond to storm events and the resulting up-stream flooding. In addition to providing a discharge-flow-depth relationship, the model will predict the required parameters to make estimates of residence time and down-stream concentration distributions using Equations (1) and (3). The validity of the frameworks assumed in Equations (1)–(3) (NRT and ADE) are also evaluated in the context of porous dams.

2. Experimental Set-up

2.1. Physical Model and Discharge Measurements

The experimental program was designed to simulate flows and the transport of pollutants introduced up-stream of a porous dam (such as a beaver dam or leaky barrier). The porous dam was modelled physically using a rectangular Perspex plate, made porous through a series of uniformly distributed circular holes, as shown in Figure 1 and Figure 2 (Case A). The dam was 440 mm high and 300 mm wide. The holes were 5 mm in diameter, 10 mm in pitch, and distributed in 29 columns and 43 rows (giving 1247 holes in total). This corresponds to an approximate overall porosity of 19%. The dam was installed on a 20 m long, 0.3 m wide, 0.47 m deep recirculating flume, 10 m from the inlet. Water was supplied to the flume from a constant head tank with a control valve at the flumes inlet to set discharge.

The dam described above, consisting of uniformly distributed holes at a constant diameter, represents the base case where the entire dam is uniformly porous (Figure 2; Case A). In addition to Case A, three further cases were investigated where sections of the dam were made impermeable to simulate sections of a real beaver dam that are effectively impermeable due to the presence of compacted mud or rocks (Figure 2; Cases B–D). For Case B, the bottom third of the dam was porous and the top two thirds were impermeable. For Case C, the middle third of the dam was porous and the bottom and top thirds were impermeable. For case D, the top third was porous and the bottom two thirds were impermeable.

Each of the four cases were undertaken with the up-stream flow-depth at pre-defined target ‘fill levels’ up-stream of the dam (See Figure 2 and Figure 3). For each fill level, the discharge was set as the maximum possible that would maintain a steady, up-stream flow-depth at the target fill level, i.e., it represents the discharge capacity of the dam at that fill level. Target fill levels were set as y = 150, 220, 290, 360, and 430 mm (Levels 1 to 5, respectively, as shown in

Table 1. Test program (fill levels are target values).

Fill Level	Measurements
(Target Depth)	Case A	Case B	Case C	Case D
1	×1 Discharge	Not run as same as Case A	Not physically possible	Not physically possible
y = 150 mm	×3 Injections
2	×1 Discharge	×1 Discharge	×1 Discharge	Not physically possible
y = 220 mm	×3 Injections	×3 Injections	×3 Injections
3	×2 Discharge	×1 Discharge	×1 Discharge	Not physically possible
y = 290 mm	×6 Injections	×3 Injections	×3 Injections
4	×1 Discharge	×1 Discharge	×1 Discharge	×1 Discharge
y = 360 mm	×3 Injections	×3 Injections	×3 Injections	×3 Injections
5	×3 Discharge	×2 Discharge	×1 Discharge	×1 Discharge
y = 430 mm	×6 Injections	×6 Injections	×3 Injections	×3 Injections

). At each fill level, the discharge was measured volumetrically, and the flow-depth was recorded. Hence, the program delivered a discharge and flow-depth at each fill level, which could then be used to build an empirical relationship to relate flow-depth to discharge, sometimes referred to as the ‘stage-discharge’ relationship, for all cases.

2.2. Dye Tracing

Pollutant transport was investigated by injecting a fluorescent dye (Rhodamine WT) 8.8 m up-stream of the dam and measuring its temporal concentration distribution both up- and down-stream of the dam. The distance between the injection location and the measurement section exceeded flow-depth × 10 to ensure the dye was cross-sectionally well mixed [61]. At least three repeat injections were conducted for each fill level. Two Turner Designs Cyclops Fluorometers were installed to measure temporal concentration profiles of the dye as it moved through the system. The first was installed at 2.98 m up-stream of the dam to measure the dye’s initial concentration distribution. A second instrument was installed 0.99 m down-stream of the dam, to measure the concentration distribution at the dam’s outlet, giving a total reach of 3.97 m. The instruments were set at mid-depth for each fill level. Figure 3 shows a long-section schematic of the test set-up.

2.3. Test Program

Table 1 provides a summary of the test series.

As shown in Table 1, in theory there are five possible test configurations for Case A (fill level 1–5), four for Case B (as fill level 1 is the same as Case 1), four for Case C (as fill level 1 is not possible) and two for Case D (as fill level 1–3 are not possible). Three repeat injections were conducted for each test and the discharge was measured once (as discharge was measured volumetrically, the measurement is already effectively time-averaged). In addition, preliminary tests were conducted for Cases A and B that were of sufficient quality to be included in the main analysis. This is why Case A fill-level 1, 3, and 5 and Case B fill-level five have additional measurements. The inclusion of the preliminary tests means that Case A has nine discharge measurements and ×8 sets of injections (rather than the minimum of five) and Case B has ×5 discharge measurements and ×5 sets of injections (rather than the minimum of four).

3. Modelling Framework

3.1. Discharge Modelling

The experimental program described in Section 2 provides discharge and flow-depth data for the four Cases A–D. It is reasonable to expect that a simple, empirical trend could be established for each of the cases that could be used to predict discharge from flow-depth or vice versa for the system in question. However, to provide a general model that can make the same prediction for a dam with any arbitrary depth of porous or impermeable sections in contact with the flow, it is necessary to establish a general relationship between discharge and the relative depth of porous or impermeable sections in contact with the flow. Each of the four cases has three possible components making up the total flow-depth; the depth below the porous section, $y_b$, the depth in contact with the porous section, $y_o$, and the depth above the porous section, $y_a$. Figure 4 provides an example of this for Case C.

To establish a general relationship where the contribution of each relative depth is accounted for, a multi linear regression analysis was utilized using MATLAB’s in-built toolbox ‘fit’. The multi-linear regression took the form:

$$Q= {\beta }_1{y}_b+{\beta }_2{y}_o+{\beta }_3{y}_a+C$$

(4)

where $Q$ is the measured discharge (m³/s);$y_b$is the flow-depth below the porous section (m); $y_o$ is the flow-depth in contact with the porous section (m); $y_a$ is the flow-depth above the porous section (m); ${\beta }_1$, ${\beta }_2$, and ${\beta }_3$ are weighting coefficients that govern the relative contribution to flow of the depths $y_b$, $y_o$, and $y_a$, respectively; and $C$ is an intercept constant.

Equation (4) is used with the measured discharge and the measured flow-depths to provide the weighting coefficients. Once the coefficients have been established, the equation can then be used to predict discharge simply on the basis of the flow-depth (broken into $y_b$, $y_o$, and $y_a$).

3.2. Concentration Profile Modelling

As discussed in Section 2.2, dye tracer (simulating a pollution event) was injected into the test facility up-stream of the dam and its concentration vs. time distribution was measured up- and down-stream. Figure 5 shows an example of the data collected.

The goal of the model developed in the paper was to be able to predict the down-stream concentration distribution of a solute solely based on the initial distribution (or assumed initial distribution) and the depths $y_b$, $y_o$, and $y_a$. The framework adopted to predict the down-stream pollutant distribution is a standard routing solution of the Advection Dispersion Equation (ADE), as presented in Equation (3). Equation (3) will allow the prediction of a down-stream concentration distribution for a measured or assumed initial distribution based on two central parameters:

• Travel time, $\overline{t}$: The difference in centroids between the up- and down-stream distribution (t₁ and t₂ in Figure 5)
• The Longitudinal Dispersion Coefficient, $D_{xx}$: A mathematical quantification of how much the solute has spread between the up- and down-stream location

Figure 5 shows an example of a prediction made using Equation (3) with the ‘best fit’ values of $\overline{t}$ and $D_{xx}$for this example, where these parameters have been optimized to give the best possible fit of the predicted down-stream distribution compared to the measured down-stream distribution. This was undertaken using a bespoke two parameter optimization code where $\overline{t}$ and $D_{xx}$ were varied until the optimal values of $\overline{t}$ and $D_{xx}$ were found, using R² as the measure of goodness of fit. Performing this optimization to the whole data set provides experimental ‘best possible’ values of $\overline{t}$ and $D_{xx}$ from the dataset within the ADE modelling framework, and forms the foundations to build a model that can predict these values solely based on the flow-depths $y_b$, $y_o$ and $y_a$.

If the system is flowing under ‘ideal’ conditions (plug-flow), the travel time $\overline{t}$ will be synonymous with the system’s Nominal Residence Time (Equation (1)). The NRT is a simple way to predict travel time for systems that can be approximated by the ideal flow assumption. Its terms are straightforward to measure in the laboratory and the calculation is trivial. The model developed in this paper will use the NRT to predict $\overline{t}$. The validity of this assumption will be fully discussed in Section 4 and Section 5.

Similar to the analysis conducted to establish a general relationship for discharge, a further multi linear regression was performed to establish a general relationship between the dimensionless Longitudinal Dispersion Coefficient and the relative dam section in contact with the flow ($y_b$, $y_o$, and $y_a$), again using MATLABs in-built toolbox ‘fit’. This analysis took the form:

$$\frac{D_{xx}}{yu}= {\alpha }_1{y}_b+{\alpha }_2{y}_o+{\alpha }_3{y}_a+c$$

(5)

where $D_{xx}$ is the Longitudinal Dispersion Coefficient (m²/s); $y$ is the total flow-depth (m); $u$is the mean velocity (m/s);$y_b$is the flow-depth below the porous section (m); $y_0$ is the flow-depth in contact with the porous section (m); $y_a$ is the flow-depth above the porous section (m); ${\alpha }_1$, ${\alpha }_2$, and ${\alpha }_3$ are weighting coefficients that govern the relative contribution of the depths $y_b$, $y_0$, and $y_a$, respectively; and $c$ is an intercept constant.

The combination of using the NRT to predict the travel time and Equation (5) for the Longitudinal Dispersion Coefficient provides the two parameters required to make a prediction of the down-stream concentration distribution within the model.

4. Results

4.1. Flow-Depth and Discharge Relationship

The relationship between the dam’s up-stream flow-depth and discharge for the four dam configurations investigated (Cases A to D) are presented in Figure 6. The recorded values of the flow-depth and the discharge for the data presented in Figure 6 are provided in

Table 2. Results for discharge vs. flow-depth for Cases A to D.

Fill Level	Case A		Case B		Case C		Case D
	y(m)	Q(l/s)	y(m)	Q(l/s)	y(m)	Q(l/s)	y(m)	Q(l/s)
1	0.145	6.271	-		-		-
	0.147	5.074
2	0.224	12.473	0.22	8.511	0.210	1.825	-
3	0.286	19.208	0.285	11.352	0.287	7.050	-
	0.290	19.414
4	0.351	24.862	0.365	13.954	0.368	10.631	0.355	1.758
5	0.422	30.445	0.430	15.947	0.420	12.403	0.426	6.379
	0.426	33.489	0.435	16.416
	0.431	34.525

. The relationship for each case is linear and is described by a standard linear regression with R² > 0.98 for all cases (see

Table 3. Linear regression trends for relationship between discharge and depth (Cases A–D).

Case	Relationship	R2
1	Q = (0.0972 y) − 0.0088	0.99
2	Q = (0.0353 y) + 0.0010	1.00
3	Q = (0.0501 y) − 0.0081	0.98
4	Q = (0.0651 y) – 0.0213	1.00

for relationships and correlations). However, each case has a unique trend, due to the dissimilar hydraulic conditions that depend on which relative sections of the dam are porous or impermeable.

The multi linear regression analysis described in Section 3.1 was performed on the dataset to establish the relative contribution of each of the dam sections.

Table 4. Discharge—Flow-depth multi linear regressing coefficients.

${\mathit{\beta }}_1$	${\mathit{\beta }}_2$	${\mathit{\beta }}_3$	$\mathit{C}$	${\mathbf{R}}^2$
0.00939	0.09396	0.03656	−0.00766	0.99

presents the weighting coefficients obtained with the multi-linear regression analysis.

Figure 7 shows a comparison between the measured and predicted discharge using the proposed model (Equation (4) and Table 4). The model shows a strong correlation to the measured discharge (R² = 0.99) and is now a general model where the same relationship is observed for all dam configurations (Cases A–D).

Equation (4) and Table 4 provide a framework that, at a given porosity, can predict the system’s discharge from a given flow-depth, split into its constituent components$y_b$, $y_o$, and $y_a$ for any arbitrary value of these depths. This model provides an engineering tool that can accurately predict the discharge of a system, at a given porosity, simply by measuring the lengths of various porous sections that are in contact with the flow. These depths are relatively simple parameters to measure in field case studies without the need for sophisticated and expensive equipment.

4.2. Quantification of Residence Time and Solute Concentration Distribution

Figure 8 shows the system’s measured travel time $\overline{t}$ (obtained from the optimization analysis described in Section 3.2) vs. NRT (Equation (1)) for all tests (values are mean of three repeats). The entire data set has a correlation between the travel time and NRT of R² = 0.98. There are a small number of runs that are not described well (only one run, Case D, fill level 4, has a percentage error >10%), which is discussed in Section 5. This high correlation across the data set suggests the system’s travel time can be approximated well by the NRT alone, and that no significant short-circuiting is taking place despite the potential dead-zones in Cases B–D. Hence, the proposed model will adopt Equation (1) to predict the system’s travel time.

The final parameter required in the model is the Longitudinal Dispersion Coefficient. Experimental values of this parameter can be obtained from the dye tests through the optimization procedure described in Section 3.2. Figure 9 shows the relationship between the dimensionless Longitudinal Dispersion Coefficient and total flow-depth for all cases. Figure 9 suggests that the dimensionless Longitudinal Dispersion Coefficient is inversely proportional to flow-depth. However, the data is scattered and there is no obvious general trend between Cases A to D.

The multi linear regression analysis described in Section 3.2 was performed on the dataset to establish the relative contribution of each of the dam sections.

Table 5. Longitudinal Dispersion Coefficient multi linear regressing coefficients.

${\mathit{\alpha }}_1$	${\mathit{\alpha }}_2$	${\mathit{\alpha }}_3$	$\mathit{c}$	${\mathbf{R}}^2$
−0.4995	−0.8397	−0.4572	0.3803	0.70

presents the weighting coefficients obtained with the multi-linear regression analysis.

Figure 10 shows the results of the analysis. The correlation between the measured and predicted Longitudinal Dispersion Coefficient across the data set was R² = 0.70. Although the correlation was lower than the analysis undertaken for the discharge and travel time prediction, the majority of the data falls within ± a factor of 2, which is generally considered reasonable for dispersion coefficient predictions [62,63,64]. Therefore, this approach was adopted to predict the Longitudinal Dispersion Coefficient within the model, and its accuracy is fully discussed in Section 5.

The results and analysis developed in the paper provide a model for the system investigated, and a framework for more general systems, that can predict discharge, travel time, and the Longitudinal Dispersion Coefficient solely using the relative depths y_b, y_o, and y_a. Based on these three predicted parameters, the model can also predict a down-stream concentration distribution using the ADE routing equation presented in Equation (3).

The final form of the model is:

$$Q= 0.00939y_b+0.09396y_o+0.03656y_a-0.00766$$

(6)

$$\overline{t} = \frac{V}{Q}$$

(7)

$$D_{xx}=uy\left(-0.4995y_b-0.8397y_o-0.4572y_a+0.3803\right)$$

(8)

$$=\frac{Q}{A}y\left(-0.4995y_b-0.8397y_o-0.4572y_a+0.3803\right)$$

$$c\left(x_2,t\right)={{\displaystyle \int }}_{\gamma -\infty }^{\infty }\frac{c\left(x_1,\gamma \right)u}{\sqrt{4\pi D_{xx}\overline{t}}}\exp \left[\frac{u^2{\left(\overline{t}-t+\gamma \right)}^2}{4D_{xx}\overline{t}}\right]\mathrm{d}\gamma$$

(9)

5. Discussion

5.1. Model Validation

The first part of the model, used to predict the system’s discharge from the three components of flow-depth, provides a general prediction of discharge that accounts for the various possible configurations of a dam with porous or impermeable sections. The high level of correlation (R² = 0.99) demonstrates that the framework developed, presented in Equation (6), can be used to accurately predict discharge for the system in question and could easily be extended as a framework for other systems.

The more complicated and less well-defined part of the model is the prediction of the down-stream concentration profile. Both the travel time and the Longitudinal Dispersion Coefficient have lower correlations to the experimental data than discharge (R² = 0.98 and R² = 0.70 respectively), and also need to be used within the ADE framework, which contains further assumptions. For this reason, some care will be taken to fully describe the accuracy of predictions made by this part of the model.

Figure 11 shows several representative examples of a predicted down-stream concentration profile made by the model (Equations (6)–(9)) (red line) compared to the experimental data (black line). The blue line shows the best possible prediction that can be made within the modelling framework of Equation (9) when travel time and the Longitudinal Dispersion Coefficient are set at the optimized values (the optimization process discussed in Section 3.2). The green line is discussed later in this section. Figure 11 presents six profiles covering a range of correlation values from good to poor. Figure 12 summarizes the correlation values for the whole data set, presenting all correlation values for the comparisons between the model and the data (clear shapes), and correlation values between the optimized values and the data (black shapes).

The correlation values for the optimized profile (with the exception of one test, Case B, level 5) all have R² > 0.8, showing that in theory, a reasonable correlation can be achieved using the modelling framework of the ADE equation (Equation (9)). The predictions made by the ADE equation with optimized values are certainly not perfect, largely due to the Gaussian transfer function assumption [60], but a correlation of R² > 0.8 is sufficient to provide a good estimate of the concentration profile that would be acceptable for most practical applications, as shown in Figure 11. This demonstrates that any lower correlation values for the actual model (using Equations (6)–(9) for $\overline{t}$ and $D_{xx}$) shows inaccuracy in the model’s ability to predict travel time and the Longitudinal Dispersion Coefficient rather than any shortcomings in the ADE framework itself for this application.

Of all the dye injections conducted in the test program, 60% of the down-stream concentration profiles were modelled with a correlation to the experimental data of R² > 0.6. This was a correlation that the authors considered an acceptable prediction of the down-stream profile (see Figure 11a–c for examples). However, 40% of the tests were modelled with an R² < 0.6. Whilst correlations R² < 0.6 still seem to provide a reasonable approximate down-stream profile, the fit is notably less satisfactory than for cases where R² > 0.6 (Figure 11d–f).

The data presented in Figure 11 provides some insight into the reason for the poor correlation for some of the model predictions. It seems that the profile predicted by the model does not always align with the data on the time-axis, leading to the low correlation in such cases. This time-axis error accounts for almost cases where R² < 0.6. This suggests that the issue with the model for the cases with low correlation is with the prediction of the travel time, rather than the Longitudinal Dispersion Coefficient. This initially seems contrary to the earlier analysis, as the travel time was modelled with a correlation to the experimental travel times of R² = 0.98 (Figure 8), whereas dispersion coefficients were modelled to the experimental data with a far lower correlation of R² = 0.70 (Figure 10).

To investigate this, a further analysis was conducted to examine the sensitivity of the model to both the travel time and dispersion coefficient parameters. An example injection was considered that has a good overall optimized correlation (R² = 0.98 when travel time and the dispersion coefficient are the optimized values). Two hypothetical cases were then considered, one where the dispersion coefficient was held at the optimized value and a systematic, increasing error was introduced to the travel time, and then the converse, where the travel time was held constant at the optimized value and a systematic, increasing error was applied to the Longitudinal Dispersion Coefficient.

The results from this analysis (Figure 13) show that the model is far more sensitive to error in the prediction of travel time than in the Longitudinal Dispersion Coefficient. In fact, a factor error of just 1.2 (i.e., $\overline{t}=\overline{t}\times 1.2$) leads to a drop in correlation from R² = 0.98 to R² < 0. For the Longitudinal Dispersion Coefficient, a factor of almost 100 (i.e., $D_{xx}=D_{xx}\times 100$) is required to produce a similar drop in correlation. The difference is so large Figure 13 had to be split into two scales, (a) and (b), to demonstrate the respective trends.

This sensitivity analysis demonstrates that the model is highly sensitive to the prediction of travel time, and that small errors in the prediction of travel time can be propagated into large errors in the final concentration profile predicted using Equation (9). This result is actually intuitive, as any error on the time-axis will shift the whole predicted profile away from the measured data.

However, despite the discovery that an almost exact prediction of travel time is required to obtain a high correlation with experimental data, this is not fatal to the model. This is because, in terms of being able to predict a down-stream concentration profile, the essential physical parameter that needs to be estimated accurately is the Longitudinal Dispersion Coefficient. The travel time is simply setting the time that the end user wants to make a prediction, whereas the Longitudinal Dispersion Coefficient will set physically how much the profile has spread at that point. Provided the prediction of the dispersion coefficient is accurate, the low correlation in some tests (due to the poor prediction of the travel time) does not necessarily mean the prediction is inaccurate, rather just that the time at which the prediction is made does not match with the occurrence of the experimental data. Provided the dispersion coefficient is accurately predicted, the profile would still be correct, just occurring at a slightly off-set time to the recorded data. For this reason, the central question is the degree to which the model is accurately predicting the Longitudinal Dispersion Coefficient.

To further investigate the accuracy of the model with respect to its prediction of the Longitudinal Dispersion Coefficient, the correlation of the model to the data when the travel time is held at the optimized value (i.e., assuming travel time is modelled perfectly) and the dispersion coefficient is predicted using the model (Equation (8)) was considered. Figure 14 presents the correlation value for the whole data set, and the green profiles in Figure 11 show example down-stream predictions using the model in this form. From Figure 14 it can be seen that, with the exception of one case, all correlation values are now R² > 0.6, and from Figure 11 it can be seen that they all make reasonable predictions of the down-stream data. This demonstrates that the model works well in terms of predicting the Longitudinal Dispersion Coefficient, which suggests the high performance of the model in predicting the down-stream concentration profile. However, a more accurate prediction of the travel time is required to improve the overall prediction of the model, given how sensitive the model is to error in travel time.

The sensitivity analysis presented in Figure 12, Figure 13 and Figure 14 shows that the main issue in poor correlation for the final concentration profile is due to small errors in the travel time prediction, combined with the model being extremely sensitive to the travel time parameter. This was initially unexpected, because this paper has demonstrated that a good prediction for the travel time can be made by estimating the system’s NRT, and that the correlation between an estimate of the NRT and the measured travel time from the dye tests was high (R2 = 0.98; see Figure 8). The NRT is also a relatively straightforward parameter to measure and calculate for a laboratory system; it only requires knowledge of the dimensions of the flume and an accurate measure of the flow-depth, which is normally trivial to measure accurately when the system is close to uniform flow conditions.

However, one of the major challenges of the experimental set-up used for this study was the ability to measure the flow-depth across the whole reach. Dye was injected 8.8 m up-stream of the dam, and the up-stream concentration profile was measured 2.98 m up-stream of the dam. Ideally, the down-stream profile would have been measured immediately down-stream of the dam. However, when conducting the tests, super-critical flow often occurred immediately down-stream of the dam, which, due to the presence of large air bubbles, did not allow for accurate measurement of the dye’s concentration. It was therefore necessary to position the down-stream concentration instrument approximately 1 m down-stream of the dam, after a hydraulic jump had occurred, so the flow was sub-critical and the concentration could be measured accurately. This meant that whilst measuring the up-stream flow-depth was trivial, the down-stream flow-depth was more challenging due to the presence of rapidly varying flow. Although it was possible to measure the flow-depth at the location of the down-stream concentration instrument, the presence of super-critical flow, a hydraulic jump and then sub-critical flow meant that the values measured were not fully representative of the actual flow-depth between the dam and the down-stream concentration instrument.

Whilst attempting to account for the uncertainty in the down-stream volume, it became clear that the most accurate residence time over the whole data set could be obtained by neglecting the down-stream volume altogether, and simply considering the volume up-stream of the dam. This worked well for the majority of cases, as the volume down-stream of the dam was often near-negligible (often far below for almost all cases 10%) compared to the volume up-stream, and the down-stream volume was so uncertain that any estimate of it led to errors. This approach led to a relatively accurate prediction of the NRT relative to the dye’s travel time (R² = 0.98). However, as discussed above, cases where the down-stream volume was not negligible, thus causing a small error in the NRT, led to large errors in the final concentration profile prediction, due to the sensitivity of the model to the travel time prediction. In future experiments, setting a control weir to maintain sub-critical flow down-stream of the dam would be considered as well as, if possible, measuring the down-stream concentration profile immediately down-stream of the dam. Alternatively, taking longitudinal high-resolution flow-depth measurements if the flow is rapidly varying would enable a more accurate prediction of the down-stream volume.

Despite the few low correlation relationships discussed above, results obtained justify that the predictions made by the model proposed in this paper provide a reasonable engineering approximation of discharge, travel time, the Longitudinal Dispersion Coefficient and the down-stream concentration profile for the majority of cases presented. The model provides a framework for predicting these parameters for porous dams with a combination of porous and impermeable sections.

5.2. Limitations and Further Work

The authors consider this work to be a pilot study and proof of concept that could easily be extended to larger scale dams at any dam porosity, or even with porosity as a further variable. The model dam is an intentionally simplistic representation of a real beaver dam, designed to simulate the essential flow parameters whilst being simple to characterize geometrically for modelling. The authors consider this study to be the first step on the ladder in the process of fully understanding the hydraulics and dispersion characteristics of such structures, which in reality are extremely complex. The modelling will also have wider applications to many other forms of porous dam. The linear regression analysis conducted to supplement the experimental work provides a simple engineering tool to make predictions for such systems on the basis of simple to measure parameters. Further work is required to develop a more sophisticated model to predict these essential parameters on a purely theoretical basis.

The modelling was undertaken in a laboratory flume, which is clearly not a full-scale model. The authors feel that it is essential to quantify these processes in a small-scale, controlled environment before the analysis is extended to larger scale models that more realistically represent real beaver dams. It is worth noting, however, that whilst the model dam is considerably narrower than a real beaver dam, the depth of the dam and pore sizes are comparable to full scale structures (at maximum flow-depth the model is of the order of 30–40% of a typical beaver dam).

The physical model of the dam used in this work is a 2D representation of a beaver dam. This simplification allows for easy characterization and helps to isolate the analysis of the effects that the front-facing area of the dam can generate on the flow. However, in reality beaver dams are 3D structures with a significant longitudinal length. To obtain full understanding of such systems, further work that considers this needs to be undertaken.

The experimental set-up in the work quantifies pollution mixing between two sites up- and down-stream of the model dam. This standard method of dye tracing between two sites provides an integrated quantification of the mixing processes between the up- and the down-stream site. The requirement of a reasonable reach length in order to deliver good quality data means that in this work, it is not possible to differentiate between mixing up-stream of the dam, at the dam location and down-stream of the dam. At one level, this is a necessary limitation of the technique and could be addressed with further work investigating the 2D/3D flow-field in the reach. However, it is also important to note that the volumes of water immediately up- and down-stream of the dam are also affected by the dam, and form as much an essential part of the mixing processes as the dam itself. As such, quantifying mixing across the reach in the immediate vicinity of the dam is essential, especially when considering parameters such as residence time, where the effect of the structure being analyzed is primarily on the volume of water up-stream of the structure.

There is a clear gap in the literature regarding quantification of the typical porosity of real beaver dams, and this posed a significant challenge when designing the model dam. Having now conducted these tests, the authors feel that the relative flows investigated were most likely higher than would be experienced for a typical beaver dam under normal operating conditions. However, the flows presented will certainly be useful at storm levels, and hence for flooding conditions which are crucial to be managed and controlled. The authors feel a similar investigation, using the modelling framework developed in this paper, needs to be conducted at a lower porosity to more accurately model real beaver dams under typical flow conditions.

Furthermore, field studies could also provide essential datasets to optimize the model suggested and help to identify areas for improvements. Recent studies [65] have demonstrated the impact that Eurasian beavers have on natural environments, and some cases are also characterized by ponds and beaver dams in series, which could have different consequences on the spread of pollutants down-stream as well as on the attenuation of flood flows. Additionally, other studies [66] also suggested that beavers bring items into the aquatic systems (e.g., corn, treated lumber) which may reduce water quality down-stream of the wetlands they create. It would be interesting to introduce these variables into the model proposed for a deeper understanding of the entire scenario.

6. Conclusions

This paper presents a novel data set that investigates the relationship between a dam with a combination of porous and impermeable sections and the dam’s fundamental hydraulics and pollutant dispersion characteristics. A semi-empirical model was developed that, for a given flow-depth (broken down into sub-depths for each porous and impermeable section), can predict the system’s discharge, residence time, Longitudinal Dispersion Coefficient, and down-stream concentration distribution for a given pollution event.

The conclusions that were drawn from this data set are that, under the flow conditions presented in this paper:

• A general relationship exists between flow-depth and discharge for dams with both porous and impermeable sections, provided a breakdown of the relative depths is known and accounted for.
• The residence time of dams with both porous and impermeable sections can be reasonably estimated using the Nominal Residence Time.
• A general relationship exists between flow-depth and the dimensionless Longitudinal Dispersion Coefficient for dams with both porous and impermeable sections, provided a breakdown of the relative depths is known and accounted for.
• The Advection Dispersion Equation can be used to predict a reasonable estimate of the down-stream concentration distribution of solutes for such dams.

The first part of the model provides a simple method/framework to predict the discharge of a stream obstructed by a porous dam, simply by measuring the relative depths of porous and impermeable sections. The second part of the model provides the ability to predict pollutant travel-time and dispersion characteristics, and allows for the prediction of how a given dam would perform in the event that pollutants were introduced at some location up-stream. This allows quantitative predictions to be made when considering the potential for such systems to aid in pollution treatment.