Modelling the Impact of Leaky Barriers with a 1D Godunov-Type Scheme for the Shallow Water Equations

: There is increasing interest in distributing small-scale interventions across the landscape as an alternative means of reducing ﬂood risk. One such intervention, the leaky barrier, is introduced in channels to slow down high ﬂows and encourage temporary storage on the ﬂoodplain. While these barriers have been implemented widely, there is still resistance to their use at the scales required to impact signiﬁcantly on ﬂood risk, at least partially due to an evidence gap. In particular, there is no standard method for representing leaky barriers in hydraulic models. This study sets out a methodology for developing mathematical models which capture the hydraulics of leaky barriers accurately, allowing key questions about their combined behaviour in catchments to be answered. A 1D Godunov-type scheme is set up and leaky barriers incorporated with internal boundary conditions. This model is tested against benchmarks from the literature and new steady-state data, and then run predictively on transient cases. The method will help to answer key questions about the optimal leakiness of small-scale interventions, the limits to their usefulness, and how combinations of barriers may or may not cause synchronisation problems when the effect of multiple barriers is


Introduction
Flooding is one of the great issues of our time, representing the world's most damaging environmental hazard [1][2][3]. The pattern of land use change in both the developed and developing world over the last century has been a key driver for increasing the risk of flooding. Most notably, urban development has had a widely recognised impact, representing a direct change from high infiltration to low infiltration [4]. Moreover, post-Green Revolution intensification of farming has increased surface runoff pathways and soil erosion, as well as introducing diffuse pollution in the form of nutrients and pesticides into waterways [5]. These impact not only on flood risk, but on numerous environmental factors including water quality, river morphology, habitats and biodiversity. While these changes are well researched and understood [6][7][8][9], solving the problems associated with them remains a major challenge. In addition, climate change poses a further threat to urban and rural catchments via wetter winters and more intense summer storms [10,11].
These current and potential future threats provide a backdrop for increasing interest in moving away from large traditional engineered structures to large numbers of small-scale interventions distributed across the landscape with a focus on approaches that mimic or work with natural estimated with 1D hydraulic models [45] and hydrological GIS tools [52]. This approach assumes the barrier is empty before the flood event, and so would not be effective at modelling multi-peak events. Also, it neglects any losses through the barrier walls, and it does not consider the location of the barriers. These problems mean that the modelled effects of the interventions are likely to be overestimates. Ideally, there would be a way to include leaky barriers in the hydrological and hydraulic models themselves.
A popular method for including leaky barriers in a model is to change parameters at the location of the barriers. For example, online barriers are often represented as increased values of Manning's n, both in hydrological models [36,53] and hydraulic models [54,55]. However, even if these parameter values come from the field, it is inappropriate to apply them to other sites or flow conditions [56]. Moreover, an increased Manning's n will impact low flows as well as high flows, and so this would not be a good model of barriers that let low flows through unimpeded. In fact, an increased Manning's n will have little effect on steeper slopes [36,57], but a leaky barrier is going to provide an obstruction either way.
A more physically-based way to model online barriers is by altering the channel geometry. For example, they have been modelled as partial blockages [47] and reduced cross-sections based on field data [58], both in 1D hydraulic models. This has a different effect to changing the roughness. Indeed, while increasing the roughness reduces velocities, decreasing the cross-sectional area increases velocities [58].
Leaky barriers can also be modelled by comparing them to well-understood hydraulic structures. For example, weir equations have been used to study the behaviour of different barrier designs-full brow, rectangular slot, V-notch, inverted V-notch, letterbox slot-as part of a hydrological model describing the ponded water behind each barrier [59]. Another option is to model leaky barriers as sluice gates with weir overflow, as was done to route the runoff from Dynamic-Topmodel [28] and as part of a simple network model [35]. While these representations have more potential to replicate the real-life behaviour of online and offline barriers, they were not validated with data. We are constrained by the theoretical and empirical equations that already exist.
There are limited detailed measurements of what happens at leaky barriers, and so there is no agreed-upon way to represent leaky barriers in models. Without further investigation, it is impossible to know if real leaky barriers will behave as high values of Manning's n, reductions in cross-sectional area, or weirs/sluice gates. To get a clearer idea of how leaky barriers behave in isolation, they have to be investigated in isolation. In particular, a recent review paper [56] has recommended that the hydraulic structure method be investigated in more detail. The aim of the present study is to develop a methodology to model leaky barriers in catchments, using a combination of physical and mathematical models. A 1D hydraulic model that includes leaky barriers as hydraulic structures is developed and then tested against experimental data obtained from a controlled hydraulic flume. The results suggest that internal boundary conditions provide a promising method for including leaky barriers as hydraulic structures in mathematical models. Figure 1a shows how the physical model of the leaky barrier consists of a single horizontal aluminium sheet with a gap underneath. Some flow will always pass under the barrier and, in high flows, there will be flow over the barrier too. In the terminology from the literature, the barrier is like an underflow jam [60]. Such barriers have been modelled by adding together the sluice gate and rectangular weir equations [28,35].

Conceptual Model
The physically-based models of [61] suggest that there are five possible operating stages for a combined gate and weir structure, as illustrated in Figure 1c. These are: 1. Gate operating freely, no weir overflow 2. Gate submerged by tailwater, no weir overflow 3. Gate and weir operating freely 4. Gate submerged by tailwater, weir operating freely 5. Gate and weir submerged by tailwater In stages 1-2, the barrier acts as a sluice gate, while in stages 3-5, the barrier acts as a combined sluice gate and rectangular weir. The differences within these two groups depend on the depth of the tailwater and the location of the associated hydraulic jump. Different discharge equations are used if a gate is operating freely, as in stages 1 and 3, or with its vena contracta submerged by a hydraulic jump, as in stages 2, 4, and 5 [62]. Similarly, different discharge equations are used if a weir is operating freely, as in stages 3 and 4, or if it is submerged, as in stage 5. Although in reality the gate and weir flows will interact [63][64][65], the interaction is ignored in the present study, and the gate and weir discharges are simply added together.

Flume Experiments
Physical modelling was undertaken in a hydraulic flume at Newcastle University. The flume has length 17.95 m, width 0.294 m, and slope 1:1600. Calibration of the flume gave a Manning's n of 0.009. A leaky barrier was placed at x = 11.14 m along the length of the flume, using slots on either side of the flume to fix it in place, as shown in Figure 1a. These slots reduced the barrier width to 0.264 m. Let the depth of the base of the leaky barrier be a 0 and the top a 1 , as in Figure 1b. Four barrier configurations were tested, each for a range of steady-state flow conditions: "Free flow" refers to stages 1 and 3, or scenarios where the leaky barrier is not touched, and "submerged" refers to stages 2, 4 and 5, where the downstream tailwater affects the flow through the barrier. Depth measurements were taken at six fixed points along the flume, two upstream of the barrier and four downstream: where the subscript denotes the distance in metres along the flume. For (c) and (d), two additional depth measurements were taken for the free-flow experiments, one either side of the downstream hydraulic jump. For each of the 55 flume experiments, the steady-state discharge was measured using a manometer (q man ). For (a), (b) and (d), a velocity meter was also used to calculate the discharge (q vel ). It was found that the manometer and velocity-meter discharges did not agree. The equation q vel = 0.00286 + 1.08 · q man (2) describes the relationship for (d) with R 2 = 0.99, as plotted in Figure 2. In case (c), there were no q vel data. However, the velocity meter was known to give more accurate results than the manometer. Therefore, for consistency, the regression Equation (2) was used to generate q vel values from measured values of q man for all cases. This transformed discharge is assumed to be an unbiased measurement for the remainder of the paper.  Figure 2. The relationship between velocity-depth and manometer discharges, with regression for (d).

Godunov-Type Scheme
The governing equations are the 1D shallow water equations, a set of nonlinear, hyperbolic partial differential equations derived from the Navier-Stokes equations with the assumption that vertical acceleration is negligible [66]. They can be written in conservative form as

∂U ∂t
where the vector of conserved variables, flux vector, and source term vector are respectively, and h is the water depth (m), u the depth-averaged velocity (m/s), g the gravitational acceleration (m/s 2 ), S f the friction slope, and S b the bottom slope. A finite volume Godunov-type scheme was chosen because it can deal with rapidly varying and transcritical flow, as well as flood waves represented as discontinuities in the water level. One splitting scheme for the inhomogeneous Equation (3) is the explicit time-marching formula where U n i is the average value of U in cell i at time t n [66]. However, the flexible semi-discrete method is used to obtain a second-order accurate scheme. This is when the space and time discretisations are considered separately [67]. Discretisation in space but not time leads to the ordinary differential equation For second-order accuracy in time, a second-order Runge-Kutta method is used to solve (6) as in [68]. where This is equivalent to (5), but second-order accuracy in time is obtained with instead. To use this Equation (9), the time step ∆t, interface fluxes F n i+1/2 and F n i−1/2 , and source term S n i need to be found. For stability, each time step ∆t is obtained using the Courant-Friedrichs-Lewy condition where a = gh is the celerity and C ∈ (0, 1] is the Courant number [66]. The value of C is fixed at the start of the simulation, although the first time step is set to be ∆t = 10 −5 s. This is because, at the start of a dam break simulation, there is no contribution to the wave speed from the particle velocity u. However, the velocity will immediately stop being zero, and so the first time step could be so large that it leads to instabilities [67].
The HLL Riemann solver of [69] is used to estimate the flux through the cell interface. Suppose the conserved variables to the left of a given interface are U L and those to the right are U R as shown in Figure 3, where O is a local origin located at the cell interface. Then the flux through the interface is approximated by where S L and S R are approximate wave speeds in the Riemann problem, estimated here using the adaptive method of [70]. In this method, an initial guess is given for the depth in the star region using a two-rarefaction approximation: If this depth h 0 is less than h L and h R , then the two-rarefaction approximation is kept to get h * = h 0 . Otherwise, a two-shock approximation is used to get where for K = L, R. Now, if h * > h L , then the left wave is a shock and otherwise Similarly, if h * > h R , then the right wave is a shock and otherwise For second-order accuracy in space, the left and right Riemann states U L and U R are found using MUSCL-type extrapolation [71] with a minmod slope limiter as in [68]. For instance, at the interface between the i th and (i + 1) th cells, the left and right Riemann states are given by with the gradient ∇U i limited by where the minmod slope limiter is defined as This requires two ghost cells at each boundary.
Only the homogeneous case was taken into account when finding the flux, even though the conservation laws (3) contain source terms. This is because the bed slope source term is unsplit from the equations, and so is incorporated into the time-marching Formula (9) at the end of each time step. A central difference approximation gives as in [68]. In contrast, the friction source term S f is split into an ordinary differential equation solved implicitly at the start of each time step with as in [68], where R is the wetted perimeter. This is stable for high values of Manning's n, which is useful when comparing different ways to model leaky barriers.

Leaky Barriers as Internal Boundary Conditions
Suppose there is a leaky barrier located exactly at the interface between the cells i and i + 1, that is, at the interface i + 1/2. The interface flux F n i+1/2 needed for the time-marching Formula (9) has two components: the mass flux hu = q and the momentum flux Unless the water is shallow enough to flow unimpeded below the barrier, these fluxes cannot be found using the Riemann solver (11)-there is a leaky barrier obstructing the flow and so the shallow water Equation (3) are no longer valid. There are two steps for solving a similar problem involving a sluice gate [72]. First, find the mass flux q through the hydraulic structure using a well-understood steady-state discharge equation. Second, find some depth h to substitute into (27) to get the momentum flux. However, the leaky barrier to be studied (Figure 1a) is not a well-understood hydraulic structure. The approach taken here is to assume that there are five possible operating stages for the barrier as outlined in Section 2.1. Let the depth of the base of the leaky barrier be a 0 and the top a 1 , and assume that the flow is from left to right (Figure 1b). First, assume that the gate and weir are operating freely. The gate flow is given by the classical equation Here, the discharge coefficient C g is a function of the upstream depth h L , where C c is the contraction coefficient, as described by [72]. The weir flow is given by where C w is a constant. Ignoring any interactions between the gate and weir, the combined flow is given by If the gate is operating freely, then there is supercritical flow downstream of the barrier, as illustrated in stages 1 and 3 of Figure 1c. The depth of this supercritical flow, denoted h sup , is calculated iteratively using an energy balance either side of the barrier, as in [73]. Following [72,74], if h R is less than the conjugate depth of h sup , then there is sufficient momentum to push the hydraulic jump downstream and the gate is operating freely. Otherwise the gate is submerged (stages 2, 4, and 5). In these cases, the submerged flow equation of [75] is used to find a reduced value of C g , where h h L = and Here, the discharge coefficient C g is a function of h R as well as h L . This is the most suitable submerged discharge equation when there are no calibration data [72,76]; however, it does add a degree of uncertainty due to its purely theoretical nature. Additionally, if h R is deeper than the critical depth over the weir, then the weir is also submerged (stage 5). In this case, the submerged flow equation of [77] is used to modify q w , multiplying it by where m = 0.185 and n = 1.5, to make q w also a function of h R as well as h L . The discharge Equation (31) provides the mass flux hu = q, which is the same leaving the left-hand cell as entering the right-hand cell. The momentum flux (27) is more difficult as it needs a value for h as well as q, and h is different in each cell due to the discontinuity created by the leaky barrier. For the momentum flux from the left-hand cell, the depth h in (27) is h L . For the flux into the right-hand cell, the depth h in (27) is h sup for stages 1 and 3, and h R for stages 2, 4, and 5.

Establishing Discharge Coefficients and Energy Losses
The internal boundary condition used in the 1D model is based on equations with empirical parameters in Supplementary Materials. These parameters were found by analysing the measured depths immediately upstream and downstream of the barrier. First, as the slots either side of the barrier reduced the barrier width to 0.264 m, the discharge Equation (31) was multiplied by the ratio 0.246/0.294 ≈ 0.84. The discharge equation requires two more coefficients: the gate contraction coefficient C c and the weir discharge coefficient C w . For each of the 52 free-flow experiments, the depth measurements h 11 were plugged into the discharge equation, with this predicted discharge then visually compared against the measured discharge as in Figure 4. The coefficients C c = 0.70 and C w = 0.70 were found to give the best fit with the measurements, and so they were used in the numerical model.  Second, in free-flow conditions, an energy balance is required to find the depth of the supercritical flow downstream of the leaky barrier, h sup . Bernoulli's equation gives Note that the change in the channel bed elevation ∆z b is not included in (36) since we are considering values at the cell interface i + 1/2. The effect of the bed slope is taken into account in the source term. If we ignore energy loss (i.e., assume ∆E = 0) and solve (36) iteratively for h sup , the water exits the leaky barrier too quickly in the 1D model, resulting in the hydraulic jump occurring further downstream than observed. However, the upstream backwater effect is captured well, indicating that the problem is not in the value of q, which determines the mass flux. This is demonstrated in Figures A1-A8 in the Appendix A. Therefore, a method for estimating ∆E in (36) was required. Incorporating energy losses into the momentum flux like this is reported elsewhere in the literature, for example, for 1D flow around islands [78].
To find a function for these energy losses, they were first approximated for each of the 52 free-flow experiments using Bernoulli's equation: where q is the measured discharge, h us is the depth measurement h 11 , and h ds is a depth measurement downstream of the leaky barrier, generally h 12 . In some cases h ds was taken to be h 11.15 or h 11.30 instead, depending on the position of the hydraulic jump and weir nappe. Note that ∆z b is included in (37) because the depth measurements were taken up to 0.3 m away from each other. This analysis resulted in the regression equation ∆E = 0.012 − 0.362 · a 0 + 0.205 · h us (38) which gave R 2 = 0.48, plotted in Figure 5. This is included in the 1D model by replacing h us with h L in (38), and then plugging ∆E into (36)

Benchmark Tests
The 1D scheme without a leaky barrier was tested to check if it correctly solved the shallow water Equation (3). These benchmarks included two steady-state and one transient case. Then the 1D scheme with a leaky barrier operating as a sluice gate (stages 1 and 2) was tested against three further benchmarks. Unless otherwise noted, a Courant number of C = 1.0 was used.

Subcritical Flow
To verify that the scheme can deal with topographical variations, the first benchmark was steady-state subcritical flow over a bump [79]. Here, the domain of length l = 25 m has its bed given by with the upstream boundary condition q(0) = 4.42 m 2 /s and downstream boundary condition h(25) = 2 m. The initial conditions are set to be q(x) = 4.42 m 2 /s and h(x) = 2 m everywhere, and Figure 6 shows that they settle out to the correct steady state. Velocity (m/s)

Transcritical Flow
To verify that the scheme can deal with transcritical flows and friction represented by Manning's n, the second benchmark was a steady-state transcritical flow case from [79]. Full derivations of the benchmark can be found in [80] and a programmed example in [81]. Figure 7 shows that, when the simulation settles to a steady state, there is only a small error in the location of the hydraulic jump.

Riemann Problem
To verify that the scheme can capture shocks, the third benchmark was a Riemann problem from [66], with discontinuous initial conditions

Sluice Gate Riemann Problems
The internal boundary condition presented in Section 2.3.2 generalises the numerical method of [72]. In the same paper [72], analytical solutions to Riemann problems where a sluice gate separates the two Riemann states were also derived. As the leaky barrier can operate as a sluice gate (stage 1 and 2), these benchmarks were used to test the 1D model. In all three cases, the gap under the barrier was set to be a 0 = 0.2 m. The initial conditions for discharge were set to be q(x) = 0 everywhere and the depth set to be where h R = 0.002, 0.2, 0.6 m for each case respectively. Figure 9 shows that the wave structure, speeds, and strengths are accurate in all three wet-bed cases, but the first test case requires a smaller Courant number to obtain accurate results. This demonstrates the shock-capturing capabilities of the scheme under transient conditions when leaky barriers are included.

Steady-State Flume Experiments
The 1D model should be able to recreate longitudinal profiles for each of the 55 steady-state flume scenarios. This was tested by using the parameters from Section 2.3.3 for the internal boundary condition, the measured discharge as an upstream boundary condition, and the depth measurement h 17.8 as a downstream boundary condition. As the grid size was set to be ∆x = 0.1, the leaky barrier was located at x = 11.1 m, four centimetres away from its real placement. The Courant number was set to be C = 1.0 and the model run until it reached a steady state.
The water surface profiles of eight representative examples are displayed in Figure 10, with full results in Figures A1-A8 in the Appendix A. First recall that the internal boundary condition introduces two new pieces of information into the domain: the discharge through the barrier and the supercritical depth downstream of the barrier. The first piece of information travels both upstream and downstream from the barrier, whereas the second piece of information travels only downstream from the barrier. Generally, the profile upstream of the barrier is simulated well by the model, reflecting the lack of scatter in the discharge equation shown in Figure 4. The profile downstream of the barrier is subject to more error, reflecting the poor R 2 value for the ∆E Equation (38) and the resulting scatter in Figure 5. Despite this scatter, the downstream profile is still captured accurately by the model in most cases. Moreover, as there were no velocity measurements across the flume, the location of the hydraulic jump is a proxy for the velocity. Its accurate placement thus suggests that the model accurately simulates the velocities.
It is interesting to note that the submerged states are recreated well, considering that the two parameters of (31) were calibrated for the free-flow conditions, and it is only the equations of [75] and [77] that were used to generalise these to submerged flow. This is important because leaky barriers are installed in series, and so the backwater effect from a downstream barrier may influence the hydraulics of an upstream barrier.

Transient Flume Experiments
Finally, to demonstrate the transient functionality of the 1D model, it was used to simulate a dam break in the flume for each configuration (a)-(d). The initial conditions were with an upstream boundary condition of q = 0.06 m 2 /s. The results, shown in Figure 11, are similar to video footage from the flume experiments. However, there is currently no quantitative data to compare them against. Nonetheless, the hydrographs in Figure 12 show that the bigger barriers provide more of an obstruction, which is to be expected.

Discussion
Despite being a popular nature-based solution, the leaky barrier is currently not well-understood. As the behaviour of individual barriers is generally unknown, inserting them into hydrological and hydraulic models involves inherent uncertainties, and so the results of these larger scale models have limitations. To address these evidence gaps, the present study tested a newly developed 1D hydraulic model with experimental flume data. The results suggest that representing leaky barriers as hydraulic structures is a valid method, provided that the barrier resembles the hydraulic structure, the energy losses are incorporated into the momentum part of the internal boundary condition, and possible submergence from downstream backwaters is considered. In particular, this investigation of the hydraulic structure method addresses a specific research gap put forward in a recent review paper on modelling leaky barriers [56]. Given the current interest discussed above in installing multiple leaky barriers across catchments for flood prevention purposes, greater understanding of potential synchronisation issues is critical. The model presented here, which provides a first approximation of a range of natural and designed barrier types including log jams, woody debris and beaver dams, has the potential to address this issue by putting a large number of leaky barriers into a 1D hydrodynamic model of a channel network.
More generally, the approach taken in this study, using a combination of physical and mathematical modelling, has great potential to answer questions which urgently need addressing, as highlighted in the introduction to this paper. Further work is currently being undertaken to take this forwards, involving further flume experiments, 1D modelling, and 2D and 3D computational fluid dynamics (CFD) simulations. For leaky barrier designs where empirical equations already exist, such as the letterbox slot design [59] which can be modelled as a large rectangular orifice with weir overflow, the approach presented here can be applied directly. However, some barrier designs do not resemble hydraulic structures with established equations (e.g., combinations of gates, weirs and orifices). Thus further research is needed to establish discharge equations. To this end, further flume experiments are being undertaken alongside CFD simulations. CFD provides a virtual hydraulic flume, allowing a greater range of barrier designs to be investigated, including under transient conditions, and will generate much greater understanding of hydraulic behaviour. The CFD models will be tested and validated using the experimental data. The results of this work will in turn be used to inform and validate 1D models for a variety of barrier designs. Ultimately these will be applied to assess the performance of networks of barriers under transient flood conditions.

Conclusions
This study set out a methodology for developing mathematical models that both capture the hydraulics of individual leaky barriers accurately and can be applied to networks of barriers across a catchment. The method was tested for a simple leaky barrier, with the 1D hydraulic modelling results comparable to those from a controlled hydraulic flume. Further work is being undertaken so that the method can be applied to more complex leaky barriers. This will answer key questions about the optimal leakiness of small-scale interventions, the limits to their usefulness, and possible synchronisation problems when the effect of multiple barriers is aggregated.