Abstract

Culverts allow roads to safely traverse small streams or drainage ditches, and their proper design is critical to ensure a safe and reliable transportation network. A correct modelization of these hydraulic structures becomes crucial in the assessment of flood footprints or discharge peak estimation in a risk evaluation plan. The question of how to include culverts comes up frequently when assembling a hydraulic model that requires the presence of as many singular elements as possible. In this work, three different culvert integrations with the surface domain are studied and compared in the context of a 2D shallow water (SW) model. All of them are based on the Federal Highway Administration (FHWA) formulation for the culvert discharge estimation but differ in complexity and in the interaction with the numerical model for surface flow, some of them as internal boundary conditions. Several steady and unsteady validation test cases are presented and the numerical results are compared with the predictions from HEC-RAS 1D and HY-8 software. The culvert area, shape and their sensitivity to the 2D computational mesh is also analyzed.

HIGHLIGHTS

  • We present a comparison among culvert formulations coupled with a 2D surface flow model.

  • The validation is done by means of a laboratory test case.

  • The influence of the computational mesh resolution is also studied.

INTRODUCTION

A culvert is a relatively short segment of conduit that is typically used to transport water underneath a roadway, railway or other type of earthen embankment. Culverts represent an essential part of drainage networks worldwide and provide an effective solution for flowing waters to cross man-made barriers. A poorly designed culvert local network could generate upstream flooding leading to serious damage (Jaeger et al. 2019a). Traditionally, culverts are designed using rainfall data from the past to estimate flow rates. Nevertheless, extreme precipitation variations increase the probability of future flood disasters due to culvert undersizing (Chang et al. 2013; Pereira et al. 2015; Jaeger et al. 2019b). Nowadays, one of the most common guidelines for culvert design is the Hydraulic Design Series Number 5 (HDS-5) – Hydraulic Design of Highway Culvert (Schall et al. 2012), where the methodology presented is based on a semiempirical estimation of the culvert discharge capacity in terms of the difference of headwater levels. This procedure has become a standardization of the manner by which culverts are designed and constructed.

In other words, a successful culvert design depends on accurately predicting the effect that these hydraulic structure will have on the surrounding region. Typically, small culverts can be expected to cause changes in the upstream water surface elevation. It is mandatory to estimate these phenomena to ensure that changes in water elevation upstream will not adversely affect the nearest areas. The hydraulic techniques to design culverts were developed several decades ago (Norman et al. 1985). These traditional culvert hydraulic design procedures are still applicable nowadays and often combined with modern computational techniques capable of dealing with complex hydraulic/hydrologic situations in real time.

The numerical simulation represents a useful tool for sizing and designing a drainage network with culverts thanks to the great power of current computers. Instead of checking culvert modifications (size, entrance type, construction material, etc.) empirically, culverts and hydraulic structures in general can be modeled within the framework of a 2D overland flow model in order to obtain a complete simulator to predict and prevent flood situations Singh et al. (2014); Fernández-Pato et al. (2016); Nan et al. (2017); Abuzied & Mansour (2018); Jaeger et al. (2019a); Echeverribar et al. (2019).

In this work, three different culvert implementations are studied and compared in the context of a 2D shallow water flow model. All of them are based on the Federal Highway Administration (FHWA) formulation for the culvert discharge estimation following the guidelines presented in Schall et al. (2012) but use different strategies to couple the culvert element and the surface flow. The first method is the most commonly used (Flumen-Group 2010; Brunner 2016; FLO-2D Software 2018) and is to assume that the culvert inlet and outlet are included in single cells of the 2D domain as volumen sink and source terms, respectively. This implies that the volume transfer between culvert and surface flow is carried out through a single surface cell. A sophistication of this model is to consider that the volume transfer is carried out through several surface cells. In both cases, the water exchange modifies only the surface water depth of the cells involved. A third possibility is to impose the culvert discharge as one of the conserved variables of the interacting surface cells in order to preserve the momentum of the culvert water flow. All these methodologies guarantee a correct water mass conservation in the whole domain.

Regarding the surface flow, the use of distributed models for hydraulic simulation provides detailed computation of the spatial variations of the hydraulic variables of interest within the domain, such as water depth and flow velocity vector. Recent work on this topic can be found in García-Navarro et al. (2019), Xia et al. (2017), Yu & Duan (2017), Bellos & Tsakiris (2016), Cea & Bladé (2015), Liang et al. (2015), Cea et al. (2010), where the applications cover a range of test cases and practical cases of different degrees of difficulty. In this work, a 2D shallow water (SW) model is used, where the equations are discretized by means of a finite volume numerical scheme (Murillo et al. 2007; Murillo & García-Navarro 2010) guaranteeing a perfect mass conservation and a correct treatment of complex hydraulic situations, as moveable wet-dry fronts. The model is capable of calculating in any kind of mesh, although in this work triangular meshes are used due to their great capacity to adapt to irregular topographies.

In recent literature, there are some examples of 2D distributed surface flow solvers which allow consideration of 1D culvert structures within the computational domain. The majority of them base the culvert flow computation in the FHWA standard guidelines but differs on how the interaction between culvert and 2D computational mesh is treated. For instance, in HEC-RAS 2D, Iber and FLO-2D models, this communication is carried out in single 2D cells at culvert inlet/outlet. Hence, the use of cells that are large enough to span at least to the bottom of the culvert is recommended (Flumen-Group 2010; Brunner 2016; FLO-2D Software 2018). In MIKE21 the computed flow through the culvert is evenly distributed along the affected cell faces (DHI-Group 2017). On the other hand, TUFLOW offers two different modelings (1D and 2D approaches) for culvert implementation (BMT-Group 2018). 1D modeling is the preferred approach when the total structure width is less than of one or two 2D cells. In this approach, momentum is not transferred into or out of the 1D element to/from the 2D domain. However, the preservation of the velocity field across the culvert is optional. In the 2D modeling approach, momentum is transferred through the structure providing far more realistic flow patterns than using a 1D element.

The temporal discretization of the equations follows an explicit scheme, which is linked to certain restrictions in the time step in order to avoid instabilities of the numerical solution. The numerical stability in unsteady shallow flows has been a matter of recent research (Murillo & García-Navarro 2010) due to the required changes of the basic numerical scheme in order to adapt it to real situations. The size of the allowable time step to guarantee stability is fundamental when applying explicit methods for hydraulic/hydrological purposes where wet/dry fronts are dynamically established.

The structure of the paper is as follows: first a brief overview of the surface flow equations and culvert modeling methods are presented in the ‘Mathematical model’ section, together with the techniques used to obtain the numerical solution. The numerical tests are presented in the Numerical tests section classified in four groups: (1) Culvert dimensions sensitivity analysis; (2) 2D Mesh sensitivity analysis; (3) Validation with laboratory data; (4) Application to a real river meander. The conclusions reached and the proposed recommendations are summarized in the Conclusions.

MATHEMATICAL MODEL

Surface flow model

In this paper, the surface flow is formulated using the 2D shallow water equations (Vreugdenhil 1994):
formula
(1)
where
formula
(2)
is the vector of conserved variables and the superscript T denotes transpose. Here h represents the water depth and and are the unit discharges, with u and v the depth averaged components of the velocity vector u along the - and -axes, respectively. The fluxes of the conserved variables can be written as
formula
(3)
where g is the acceleration due to gravity.
We now describe the two terms on the right-hand side of (1). The term S corresponds to friction and bed slope and is defined as
formula
(4)
where are the friction slopes in the - and -directions, respectively. They are written in terms of the Manning's roughness coefficient :
formula
(5)
The term H in (1) is defined by
formula
(6)
where refers to a discharge per unit area of the cells connecting the culverts and the 2D domain.
Assuming dominant advection, (1) can be classified as a hyperbolic system. Its solution is approximated with the well-balanced explicit, first-order, upwind finite volume scheme described in Murillo & García-Navarro (2010) by integrating (1) at each cell of the computational mesh (Figure 1):
formula
(7)
where represents the fluxes in the normal direction to C. Each of the cell variables is updated by means of the Jacobian matrix of these fluxes, the eigenvalues and the eigenvectors , being the wave index. Hence, the expression for each cell in a triangular mesh can be reordered as follows:
formula
(8)
where is the cell area, the length of each cell edge and stands for the compact representation of fluxes and sources contributions at each cell edge w (Murillo & García-Navarro 2010). The minus superscript denotes the ingoing contributions arriving to each cell, according to the upwind methodology. In order to take into account that there is a finite number of p locations in the 2D domain where the culvert participates, the term within the mass source term can be formulated as follows:
formula
(9)
where if cell i corresponds to a connection cell with any culvert inlet or outlet and it is nil otherwise. represents the culvert discharge.
Figure 1

Cell connectivity sketch in a triangular unstructured mesh.

Figure 1

Cell connectivity sketch in a triangular unstructured mesh.

As the temporal discretization in (8) is explicit, the time step choice must guarantee the stability of the numerical solution. Here, is dynamically chosen as follows:
formula
(10)
where CFL is the Courant-Friedrichs-Levy number, restricted to in rectangular grids and to in unstructured triangular meshes, and
formula
(11)

The wet/dry fronts are well tracked providing stable solutions due to the use of a dynamical control of the time step size with a numerical mass error of the order of the machine precision. The use of a distributed surface flow model allows calculation of all the hydraulic variables, such as the water depth h or the flow velocities in every cell of the computational mesh.

Culvert flow computation

The procedure followed in this work for the culvert discharge estimation has been developed by the FHWA. A first computation of the culvert discharge is carried out, assuming that the culvert barrel is capable of conveying more flow than the inlet will accept (inlet control). Then, the discharge is computed again, assuming that the culvert barrel is not capable of conveying as much flow as the inlet opening will accept (outlet control). Both values are compared and the lower of the two is selected as the basis of the culvert design.

The complexity associated with inlet control, when combined with the large number of different shapes, sizes, and entrance conditions available for culverts, makes it nearly impossible to develop a single formula capable of describing the hydraulic behavior of culverts operating under inlet control. As a result, empirical coefficients are typically used to evaluate inlet control. Following Norman et al. (1985), Froehlich (2003), and Schall et al. (2012), the culvert discharge is computed in general as:
formula
(12)
where is the number of identical barrels, is a discharge coefficient that depends on the flow control and culvert geometric characteristics and is the culvert area at full section.
For inlet control calculation, is defined as
formula
(13)
being the water surface elevation at the culvert inlet, the inlet invert elevation (see Figure 2) and the discharge coefficient :
formula
(14)
where is the culvert diameter for circular culverts and the base dimension for box culverts, for mitered inlets and for all other inlets and Y, K, M, are inlet control coefficients.
Figure 2

Culvert variables description.

Figure 2

Culvert variables description.

For outlet control, is defined as
formula
(15)
being the water elevation downstream (tailwater) (see Figure 2) and the following formula is used to determine :
formula
(16)
where is the culvert hydraulic radius, is the entrance loss coefficient, is the Manning's roughness coefficient within the conduit and is the culvert length.

The culvert computation algorithm is presented within the flowchart shown in Figure 3. Once the culvert discharge is properly estimated according to the flow conditions, the surface water depth of the interacting cells are updated. From this basis, three different strategies to include culverts in a 2D model are compared.

Figure 3

Culvert computation algorithm.

Figure 3

Culvert computation algorithm.

Culvert model 1: single-cell interaction

The culvert is assumed to connect two single cells of the surface mesh, regardless of the conduit transversal dimensions (see Figure 4, upper). The computed discharge (Equation (12)) is transformed into a volume transfer , that is substracted from the inlet element and added to the outlet element assuming instantaneous water volume transmission (Brunner 2016):
formula
(17)
formula
(18)
being and the subindices referring to headwater and tailwater cells, respectively.
Figure 4

Scheme of the surface cells involved according to the culvert model. Two possible situations are depicted for model 1, one with a width similar to the cell size and another with a much wider.

Figure 4

Scheme of the surface cells involved according to the culvert model. Two possible situations are depicted for model 1, one with a width similar to the cell size and another with a much wider.

Figure 5

Culvert discharge splitting for model 3.

Figure 5

Culvert discharge splitting for model 3.

Culvert model 2: distributed interaction via the water level

In this case, the number of surface interacting cells is chosen depending on the culvert width (see Figure 4, lower). The culvert discharge is computed as in (12) by replacing the inlet/outlet water levels () by the average values in all inlet/outlet cells ():
formula
(19)
formula
(20)
where is the number of surface cells participating in the volume exchange at each side of the culvert. The volume transfer between every interacting cell i at both culvert sides is carried out as in the previous model. There is no need to force mesh breaklines or the number of cells at the inlet to match the outlet because this formulation of the culverts includes the presence as volume source terms and not as boundary conditions:
formula
(21)
formula
(22)

Culvert model 3: distributed interaction via the water discharge

The culvert discharge is computed as in the previous model and distributed among the surface cells participating at inlet and outlet (Figure 4, lower). In this case, they belong to a mesh breakline. The main difference lies in the fact that now the culvert discharge (per unit length)
formula
(23)
is split in x and y directions:
formula
(24)
being the angle between the culvert longitudinal edge and the x axis (Figure 5) and directly imposed as internal boundary conditions at the participating surface cells:
formula
(25)
formula
(26)
The surface water depth is updated by the surface flow computation.

NUMERICAL TESTS

Test 1: culvert dimensions sensitivity analysis

The case setup consists of a domain with no slope. Manning roughness coefficient is set to for all the domain and for the culvert conduit. The cross-section of the conduit is assumed to be rectangular and the invert elevations of both sides of the culvert has been set to . A constant initial water level of with zero velocity is assumed throughout the domain. The inlet discharge is set to a constant value of at the left boundary whereas a constant water level of is imposed downstream. Figure 6 shows the case geometry and the initial and boundary conditions.

Figure 6

Test 1. Geometry and initial and boundary conditions.

Figure 6

Test 1. Geometry and initial and boundary conditions.

In order to test the sensitivity of the model to changes in span () and rise () of the culvert, a total of 10 culvert cross-sections have been considered (see Table 1 and Figure 7). The same unstructured 10,046 cell mesh has been used for all the cases. Cases 1.1, 1.6–1.10 have the same cross-sectional area, as well as Cases 1.4 and 1.5 in order to evaluate separately the influence of the span and rise, keeping the same culvert area. All the cases are computed using culvert model 1.

Table 1

Test 1. Culvert span and rise values

Setup () () ()
1.1 1.0 1.2 1.2 
1.2 2.0 2.4 4.8 
1.3 0.5 0.6 0.3 
1.4 1.0 3.6 3.6 
1.5 3.0 1.2 3.6 
1.6 3.0 0.4 1.2 
1.7 1.2 1.0 1.2 
1.8 0.4 3.0 1.2 
1.9 0.2 6.0 1.2 
1.10 0.12 10.0 1.2 
Setup () () ()
1.1 1.0 1.2 1.2 
1.2 2.0 2.4 4.8 
1.3 0.5 0.6 0.3 
1.4 1.0 3.6 3.6 
1.5 3.0 1.2 3.6 
1.6 3.0 0.4 1.2 
1.7 1.2 1.0 1.2 
1.8 0.4 3.0 1.2 
1.9 0.2 6.0 1.2 
1.10 0.12 10.0 1.2 
Table 2

Test 2. Upstream water level values for all the meshes considered and relative differences with HY-8 predicted value, taken as reference

Cells
153 1.05 2.01 3.83% 
1,690 2.80 2.22 6.08% 
7,601 5.95 2.40 14.83% 
10,046 9.45 2.78 32.78% 
18,349 11.20 2.92 39.71% 
40,796 18.20 3.42 63.64% 
50,433 21.35 3.71 77.46% 
Cells
153 1.05 2.01 3.83% 
1,690 2.80 2.22 6.08% 
7,601 5.95 2.40 14.83% 
10,046 9.45 2.78 32.78% 
18,349 11.20 2.92 39.71% 
40,796 18.20 3.42 63.64% 
50,433 21.35 3.71 77.46% 
Figure 7

Test 1. Culvert span and rise comparison.

Figure 7

Test 1. Culvert span and rise comparison.

Figure 8 shows the temporal evolution of water discharge inside the culvert () and water surface elevation values at the inlet (WSEL1) and outlet (WSEL2), as well as the variations on the type of control that governs the culvert flow (inlet control or outlet control). As expected, the flow control varies more frequently when the culvert is not totally submerged, as in Case 1.2. An abrupt change in is observed in Case 1.3 at due to the dam overtopping of the upstream water level (). These hydraulic elements generate sudden changes in the bottom level when overtopping occurs, which becomes a very complex numerical situation involving the combination of both surface and culvert flows.

Figure 8

Cases 1. 1 to 1.3. Culvert discharge and WSELs.

Figure 8

Cases 1. 1 to 1.3. Culvert discharge and WSELs.

Figure 9 shows the comparison of the culvert discharge, headwater and tailwater WSEL for Cases 1.1, 1.6 and 1.7. The control type is also shown by means of a switch variable (0 = outlet control, 1 = inlet control). Most of the differences are observed during the transient phase, specially in the control type variable. Nevertheless, the same steady state is reached in all cases with minor differences in WSEL1 for the setup 1.6 due to its large rise value, which leads to unsubmerged flow conditions. The same comparison is presented in Figure 10 for Cases 1.4 and 1.5, with the peculiarity that the cross-section area is the same, rotated 90 degrees. In this case, both setups achieve exactly the same steady state for discharge and water surface elevations, but we observe notable differences in the control type variable.

Figure 9

Cases 1, 6 and 7 comparison. Full time range (upper) and detail of the inlet/outlet control changes during transient state.

Figure 9

Cases 1, 6 and 7 comparison. Full time range (upper) and detail of the inlet/outlet control changes during transient state.

Figure 10

Cases 4 and 5 comparison. Full time range (upper) and detail of the inlet/outlet control changes during transient state.

Figure 10

Cases 4 and 5 comparison. Full time range (upper) and detail of the inlet/outlet control changes during transient state.

Finally, Figure 11 shows the WSEL longitudinal profile of the full domain at steady state for Cases 1.1 and 1.6 to 1.10, which share the cross-sectional area value. The WSEL value strongly depends on the culvert width for a fixed computational mesh. The wider the culvert in relation to the cell size, the higher the upstream WSEL value.

Figure 11

Longitudinal profile of the channel showing WSEL at steady state for setups 1.1, 1.6, 1.7, 1.8, 1.9 and 1.10.

Figure 11

Longitudinal profile of the channel showing WSEL at steady state for setups 1.1, 1.6, 1.7, 1.8, 1.9 and 1.10.

Test 2: mesh sensitivity analysis

In this section, the same geometry, initial conditions and boundary conditions are used as in Test 1. However, the culvert size is fixed to and the computational mesh is modified in order to check the sensitivity to this feature of the model (Table 2). Figure 12 shows a detail of three of the meshes considered in this section (10,046, 1,690 and 50,433 cells) and the comparison with the culvert width. In order to characterize properly the relation between the culvert width and the cell size, the dimensionless parameter d is defined as follows:
formula
(27)
being the average face length of the computational mesh.
Figure 12

Test 2. Mesh comparison: 10,046 cells (upper), 1,690 cells (center) and 50,433 cells (lower).

Figure 12

Test 2. Mesh comparison: 10,046 cells (upper), 1,690 cells (center) and 50,433 cells (lower).

Results for culvert model 1

Figure 13 overlaps the culvert discharge for the three chosen meshes. During the transient state, the culvert flow is influenced by the size of the mesh compared to the size of the culvert but the same steady discharge is finally reached in all cases.

Figure 13

Test 2. Discharge evolution for the three chosen meshes (culvert model 1).

Figure 13

Test 2. Discharge evolution for the three chosen meshes (culvert model 1).

Figure 14 shows the longitudinal WSEL profile along the channel at steady state for the three studied meshes. The comparison with the results generated by HEC-RAS 1D and HY-8 Culvert Hydraulic Analysis Program, provided by the FHWA, is also plotted. We observe that upstream steady water level is strongly influenced by the value of d parameter. The coarser the mesh, the more similar the results are to those produced by HEC-RAS or HY-8 1D models.

Figure 14

Test 2. Longitudinal WSEL profile at steady state for the three studied meshes. Comparisons with HEC-RAS 1D and HY-8 are also shown (culvert model 1).

Figure 14

Test 2. Longitudinal WSEL profile at steady state for the three studied meshes. Comparisons with HEC-RAS 1D and HY-8 are also shown (culvert model 1).

In order to study more in detail the dependence of the numerical results with d, some additional meshes are added to the comparative benchmark. Table 3 shows the 7 meshes considered in this section, together with their d value, the numerical upstream steady water level and the relative difference of this value with the one provided by the HY-8 software, assumed as the reference value for the comparisons (). When plotting these values against d parameter, we see a clear linear trend (Figure 15), leading to the conclusion that large values of d generate unacceptable results for the water level upstream the culvert location. The fact that the numerical results show such dependence with the mesh resolution is totally unacceptable for a proper modeling of culvert within the framework of a 2D surface flow simulator. Hence, an improvement of this culvert model is mandatory.

Table 3

Test 2. Upstream water level values for the three culvert models and relative differences with HY-8 predicted value, taken as reference

Cells
1,690 2.80 2.217(6.08%) 2.037(2.54%) 2.036(2.58%) 
10,046 9.45 2.775(32.78%) 2.078(0.57%) 2.078(0.57%) 
50,433 21.35 3.709(77.46%) 2.098(0.38%) 2.108(0.86%) 
Cells
1,690 2.80 2.217(6.08%) 2.037(2.54%) 2.036(2.58%) 
10,046 9.45 2.775(32.78%) 2.078(0.57%) 2.078(0.57%) 
50,433 21.35 3.709(77.46%) 2.098(0.38%) 2.108(0.86%) 
Figure 15

Test 2. Linear fit to the relative error () vs. dimensionaless culvert width .

Figure 15

Test 2. Linear fit to the relative error () vs. dimensionaless culvert width .

Results for culvert model 2

In this section, results for culvert model 2 considering the 1,690, 10,046 and 50,433 cell meshes (Figure 12) are presented. Figure 16 shows the culvert discharge whereas Figure 17 shows the WSEL longitudinal profile at steady state. Unlike the previous culvert model, all the discretizations converge to the same upstream WSEL value, which matches the results from HY-8 and HEC-RAS 1D models. The dependence of the numerical results with d parameter is removed when distributing the volume transfer in several cells, showing a notable improvement in the numerical results.

Figure 16

Test 2. Discharge evolution for the three chosen meshes (culvert model 2).

Figure 16

Test 2. Discharge evolution for the three chosen meshes (culvert model 2).

Figure 17

Test 2. Longitudinal WSEL profile at steady state for the three studied meshes. Comparisons with HEC-RAS 1D and HY-8 are also shown (culvert model 2).

Figure 17

Test 2. Longitudinal WSEL profile at steady state for the three studied meshes. Comparisons with HEC-RAS 1D and HY-8 are also shown (culvert model 2).

Results for culvert model 3

Results obtained when applying the culvert model 3 to this steady case are similar to those provided in the previous sections. Both the culvert discharge and WSEL steady values are correctly reached regardless of the computational mesh (see Figures 18 and 19).

Figure 18

Test 2. Discharge evolution for the three chosen meshes (culvert model 3).

Figure 18

Test 2. Discharge evolution for the three chosen meshes (culvert model 3).

Figure 19

Test 2. Longitudinal WSEL profile at steady state for the three studied meshes. Comparisons with HEC-RAS 1D and HY-8 are also shown (culvert model 3).

Figure 19

Test 2. Longitudinal WSEL profile at steady state for the three studied meshes. Comparisons with HEC-RAS 1D and HY-8 are also shown (culvert model 3).

Results comparison

Table 3 shows the comparison in WSEL among all the culvert models considered in this work. Culvert models 2 and 3 present admissible differences () with HY-8 software, whereas culvert model 1 provides unacceptable WSEL values when using high resolution computational meshes ().

In order to emphasize the differences among culvert models from a bidimensional perspective, Figure 20 shows a detail of the culvert inlet at t = 5s for the 9 combinations of culvert model/mesh. In the light of this distributed view of the numerical results, it is clear that the more refined the computational mesh is, the more adequate are culvert models 2 and 3 with respect to culvert model 1, where the flow pattern is not realistic in 10,046 and 50,433 cell meshes.

Figure 20

Comparison among water levels provided by all the culvert models for the three considered meshes at .

Figure 20

Comparison among water levels provided by all the culvert models for the three considered meshes at .

Test 3: validation with laboratory measured data

This test aims to validate the proposed strategies using a set of experimental data specially designed for culvert model validation. The experiments were carried out in the Laboratory of Hydraulics of the University of Zaragoza in a recirculating, free-surface open flume made of methacrylate. The methacrylate flume () was long with a constant rectangular cross-section and flat longitudinal slope. A rectangular obstacle long, width and height was placed downstream the flume inlet, with a rectangular culvert in the center. The inflow to the flume was from an upstream reservoir, long and width, which was fed using a recirculation system from a recovery tank placed at the end of the channel. The inflow discharge at the upstream reservoir was controlled using a flowmeter inserted in the recirculation conduit. Two steady flow conditions were studied at and . When setting steady discharge the flow at the culvert was in free surface regime, whereas the discharge led to overtopping flow.

The water surface elevation at the culvert region was measured using a pioneering technique which allows 2D transient measurements of the flow free surface to be obtained (Martínez-Aranda et al. 2018). A RGB-D sensor (Microsoft Kinect, released in 2010) was suspended above the flume floor, ensuring a good compromise between field-of-view, 2D spatial resolution (millimeters-per-pixel) and depth-accuracy (Caviedes-Voullième et al. 2014), and covering the whole culvert region. A sketch of the experimental setup is shown in Figure 21. This device provides a sequence of RGB + Depth VGA binary images of the objects placed into its field-of-view. Briefly, the Kinect sensor works using the standard structured light (SL) principle, i.e. projecting an infrared pattern by means of a NIR laser diode at wavelength onto the objects. The apparent pattern deformation due to the position and shape of the objects is recorded by a monochrome NIR camera observing from a slightly different angle. This apparent deformation allows the device to produce – in hardware – a depth map for the VGA image. The 2D spatial resolution of the VGA images was and the depth accuracy was approximately for the distance between the sensor and the water surface considered in these experiments (Khoshelham & Elberink 2012). In order to allow the sensor to observe the water surface, the projected infrared pattern needs to be reflected by the flow free surface. A simple solution is to tint the water until it is quasi-opaque. In the present work, water was tinted with titanium dioxide () at a concentration of in mass. The was previously mixed with the upstream reservoir fluid and the tinted water was recirculated continuously avoiding deposition.

Figure 21

Test 3. Experimental setup for the water free surface measurements.

Figure 21

Test 3. Experimental setup for the water free surface measurements.

The device streams both the RGB and depth binary images, which are directly recorded in a solid-state disk using an ad-hoc C ++ code based on the open-source libfreenect library (OpenKinect 2012) with an acquisition rate of approximately. Each captured RGB and depth binary image was timestamped with millisecond resolution. For each experiment, approximately 90 2D binary images of the flow free surface were recorded. Once an experiment was carried out, a post-processing ad-hoc C ++ code was used to combine the stored RGB-Depth binary images into an unstructured 3D point-cloud for each measurement. The 3D point-clouds were projected onto a 2D raster grid of spatial resolution and averaged in time, in order to obtain unique uniformly distributed experimental data map of the flow free-surface position. It is worth noting that this technique allows the free surface position upstream, downstream and over the culvert obstacle to be obtained, but not inside the culvert.

For the numerical simulation, three computational meshes have been considered (2,529 cells, 11,291 cells and 149,838 cells) in order to test the sensitivity of the three culvert models to the mesh spatial resolution. Figures 22 and 23 show the measured data and the numerical results for the water depth at the steady state, for and , respectively. The results corresponding to culvert model 3 are not included as it provides the same numerical results as culvert model 2 for this test case.

Figure 22

Test 3. Validation with laboratory measured data: water depths at steady state ().

Figure 22

Test 3. Validation with laboratory measured data: water depths at steady state ().

Figure 23

Test 3. Validation with laboratory measured data: water depths at steady state ().

Figure 23

Test 3. Validation with laboratory measured data: water depths at steady state ().

Figure 24 shows the longitudinal profiles at steady state for the highest discharge () where overtopping occurs. It is worth remarking that, in this situation, the flow downstream the culvert is the sum of the flow within and above the culvert. These figures of steady profiles show that culvert model 2 produces numerical results much closer to experimental measurements than culvert model 1. Additionally, the mesh dependency is significantly reduced with culvert model 2.

Figure 24

Test 3. Validation with laboratory measured data: longitudinal profiles at steady state ().

Figure 24

Test 3. Validation with laboratory measured data: longitudinal profiles at steady state ().

As this test case has a strongly marked 2D structure, in order to quantify the quality of the numerical results for both culvert models, the differences with respect to the measured data are computed point by point along the whole 2D domain. Then, these differences are synthesized as the root mean square error (RMSE) for each culvert model in all the considered meshes (see Tables 4 and 5). In general, culvert model 2 produces significantly more accurate results, achieving 122% improvements over culvert model 1. Another conclusion is the fact that the improvements are more evident as the computational mesh becomes finer.

Table 4

Test 3. RMSE error values and % of improvement of culvert model 2 over culvert model 1 () for

Cells
2,529 3.5 16.1 8.1 98.4% 
11,291 7.5 24.7 12.3 101.0% 
149,838 15 29.5 13.9 111.4% 
Cells
2,529 3.5 16.1 8.1 98.4% 
11,291 7.5 24.7 12.3 101.0% 
149,838 15 29.5 13.9 111.4% 
Table 5

Test 3. RMSE error values and % of improvement of culvert model 2 over culvert model 1 () for

Cells
2,529 3.5 11.3 10.3 9.5% 
11,291 7.5 16.1 7.8 107.3% 
149,838 15 20.6 9.3 122.0% 
Cells
2,529 3.5 11.3 10.3 9.5% 
11,291 7.5 16.1 7.8 107.3% 
149,838 15 20.6 9.3 122.0% 

Test 4: application to a river meander

In this test case, a single meander of the Ebro river (Zaragoza, Spain) is considered. It is provided with a floodbank. The actual Digital Terrain Model (DTM) of the reach (Figure 25) and the discharge peak values have been taken from IGN (Instituto Geográfico Nacional) and CHE (Confederación Hidrográfica del Ebro) databases. The culvert system has been assumed to act as flow spillways in flood situations in order to partially absorb the river discharge peaks. Figure 25 also shows the domain delimitation as well as floodbank and culvert locations. Inflow and outflow open boundaries are also depicted. Even thought this is a synthetic case without field data, the performance of the three culvert models will be evaluated.

Figure 25

Test 4. Digital Terrain Model of the considered meander and computational domain delimitation. The floodbank and culvert locations are also shown, as well as the inlet and outlet open boundaries.

Figure 25

Test 4. Digital Terrain Model of the considered meander and computational domain delimitation. The floodbank and culvert locations are also shown, as well as the inlet and outlet open boundaries.

The domain discretization is done by means of an unstructured triangular mesh (Figure 26), locally refined with 3 different cell sizes for each domain element (see Table 6). Figure 27 shows a detail of the three refinement levels used in the computation: 54,000 cells (upper), 159,000 cells (center) and 227,000 cells (lower). Figure 28 shows the inlet discharge function.

Table 6

Test 4. Mesh construction parameters for the chosen refinements

MeshCellsOutlineInnerOuterRiverbankFloodbank
227,000 50 m 35 m 35 m 10 m 2 m 
54,000 50 m 50 m 50 m 15 m 5 m 
159,000 75 m 75 m 75 m 25 m 2 m 
MeshCellsOutlineInnerOuterRiverbankFloodbank
227,000 50 m 35 m 35 m 10 m 2 m 
54,000 50 m 50 m 50 m 15 m 5 m 
159,000 75 m 75 m 75 m 25 m 2 m 

The length in m express the maximum value of the cell face for each region of the mesh.

Figure 26

Test 4. Locally refined computational mesh.

Figure 26

Test 4. Locally refined computational mesh.

Figure 27

Test 4. Detail of the three refinement levels used in the computation: 54,000 cells (upper), 159,000 cells (center) and 227,000 cells (lower).

Figure 27

Test 4. Detail of the three refinement levels used in the computation: 54,000 cells (upper), 159,000 cells (center) and 227,000 cells (lower).

Figure 28

Test 4. Inlet hydrograph.

Figure 28

Test 4. Inlet hydrograph.

Results for culvert model 1

In this section and in the following, the numerical results obtained with the three culvert models are presented. In order to provide a visual perspective of the water depth evolution along the whole domain, Figure 29 shows the water depth values for several instants of the simulation using the culvert model 1. For a detailed comparison among culvert models and meshes, we will focus on the culvert closest to the inflow. Figure 30 shows the floodbank sectional representation and WSEL values at culvert location at , and . Note that the bed level representation slightly differs depending on mesh refinement. It is worth mentioning that all the meshes provide similar numerical results for WSEL before the culvert inlet elevation is reached (Figure 30 (upper)). However, large differences (up to ) arise at when the water starts flowing through the culvert conduit. Hence, as in the previous test cases, there exists a mesh dependence of the WSEL values with the computational mesh when using culverts model 1.

Figure 29

Test 4. Water depth evolution for , , , , and .

Figure 29

Test 4. Water depth evolution for , , , , and .

Figure 30

Test 4. Floodbank sectional representation and WSEL values at culvert 1 location at , and (culvert model 1).

Figure 30

Test 4. Floodbank sectional representation and WSEL values at culvert 1 location at , and (culvert model 1).

Results for culvert model 2

When repeating this test case with culvert models 2 and 3, the differences in water levels among meshes are reduced noticeably and can be attributed to the different bed level discretization (see Figure 31). Results for culvert model 3 are not presented due to the great similarity with culvert model 2.

Figure 31

Test 4. Floodbank sectional representation and WSEL values at culvert 1 location at , and (culvert model 2).

Figure 31

Test 4. Floodbank sectional representation and WSEL values at culvert 1 location at , and (culvert model 2).

Result comparison

In order to perform a detailed comparison among numerical results, Figure 32 shows the water level differences between mesh pairs for culvert model 1 and 2. Culvert model 2 clearly reduces the mesh dependence for all the comparisons performed. A quantification of these differences is presented in Table 7, that summarizes the average differences in WSEL and flow velocity modulus U along the floodbank longitudinal profile for culvert models 1 and 2. Diff. (%) stands for the relative difference between both culvert models. A maximum improvement of in WSEL and in U is achieved when using the culvert model 2.

Table 7

Test 4. Average differences in WSEL along the floodbank longitudinal profile for culvert models 1 and 2

Culvert model 1Culvert model 2Diff. (%)
 0.08930 0.04337 206% 
 0.12300 0.04517 272% 
 0.09677 0.05907 164% 
 0.64289 0.42566 151% 
 0.85957 0.69368 124% 
 0.55455 0.46988 118% 
Culvert model 1Culvert model 2Diff. (%)
 0.08930 0.04337 206% 
 0.12300 0.04517 272% 
 0.09677 0.05907 164% 
 0.64289 0.42566 151% 
 0.85957 0.69368 124% 
 0.55455 0.46988 118% 

Diff. (%) stands for the relative difference between both culvert models.

Figure 32

Test 4. Water level differences between mesh pairs for culvert models 1 and 2.

Figure 32

Test 4. Water level differences between mesh pairs for culvert models 1 and 2.

CONCLUSIONS

A proper representation of culverts in 2D surface flow models is essential when dealing with flood event simulations, since those hydraulic structures artificially convey the water flow and modify the runoff capacity of natural and urban catchments. This paper presents three different culvert numerical implementations within the framework of a 2D distributed shallow water flow simulator. Several test cases have been presented in order to analyze the influence of the relative culvert-mesh size on discharge and water surface elevations computation, paying special attention to the upstream values. In the light of the numerical results, the conclusions of this paper can be summarized in the following points:

  • The culvert relative dimension to the grid size sensitivity analysis provided some insights through a set of ten test cases with different culvert cross-sections. In general, for a fixed computational mesh, the upstream water level value strongly depends on the culvert width when using the culvert model 1. The wider the culvert in relation to the cell size, the higher the upstream WSEL value. On the other hand, the setups with common cross-sectional areas and flow regime achieve exactly the same steady state for discharge but present notable differences in the control type variable during the transient state.

  • Regarding the mesh sensitivity analysis with a fixed culvert cross-section, we observe that the upstream steady water level is totally influenced by the mesh size compared to the culvert width when using the culvert model 1, as in the previous test cases. The coarser the mesh, the more similar the results are to those produced by HEC-RAS 1D or HY-8. Nevertheless, the steady discharge is correctly matched. A total of 7 meshes are used in order to develop some insights on how the ratio between culvert width and cell size ( parameter) influences the upstream water level, achieving a linear fit between the simulation error and d values. When performing this mesh sensitivity tests using culvert models 2 and 3, the mesh dependence is totally avoided, providing consistent numerical results. In addition, the logical 2D flow patterns are preserved regardless of the mesh resolution or culvert dimensions.

  • A laboratory test has been performed in order to validate the proposed strategies for culvert implementation. The experiments were carried out in a recirculating, free-surface open flume with a rectangular obstacle provided with a rectangular culvert. Two steady flow conditions were studied, leading to a free surface regime and overtopping flow, respectively. The acquisition of the water depth measurements was carried out by a RGB-sensor. For the numerical simulation, three computational meshes have been considered (2,529 cells, 11,291 cells and 149,838 cells) in order to test the sensitivity of culvert models 1 and 2 to the mesh spatial resolution. In general, culvert model 2 produces significantly more accurate results, achieving 122% improvements over culvert model 1. Another conclusion is the fact that the improvements are more evident as the computational mesh becomes finer.

  • The application to a real topography has been presented in the third test case. A river meander provided with a floodbank and several culverts has been considered. Again, a dependence of the upstream water level with the computational mesh arises when using the culvert model 1 reaching differences among meshes of larger than the ones obtained with culvert models 2 and 3. The diferences among meshes in flow velocity models are also reduced by up to .

In general terms, the use of a dynamic choice of interacting cells between culvert and surface domain provides accurate numerical results, consistent with the water surface elevations predicted by 1D models.

ACKNOWLEDGEMENTS

The present work was partially funded by the Aragón Government through the Fondo Social Europeo. It has also been supported by the Research Project PGC2018-094341-B-I00 funded by the Spanish Ministry of Science, Innovation and Universities (MICINN). Additionally, Mario Morales-Hernández was partially supported by the U.S. Air Force Numerical Weather Modeling Program.

DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

REFERENCES

REFERENCES
BMT-Group
2018
TUFLOW Classic/HPC User Manual, Build 2018-03-AD
.
Brunner
G.
2016
HEC-RAS River Analysis System, 2D Modeling User's Manual Version 5.0
.
U.S. Army Corps of Engineers
,
Davis, CA
,
USA
.
Chang
H. K.
Tan
Y. C.
Lai
J. S.
Pan
T. Y.
Liu
T. M.
Tung
C. P.
2013
Improvement of a drainage system for flood management with assessment of the potential effects of climate change
.
Hydrological Sciences Journal
58
,
1581
1597
.
DHI-Group
2017
Mike21 Flow Model FM, Hydrodynamic Module User Guide
.
Echeverribar
I.
Morales-Hernández
M.
Brufau
P.
García-Navarro
P.
2019
Use of internal boundary conditions for levees representation: application to river flood management
.
Environmental Fluid Mechanics
19
,
1253
1271
.
Fernández-Pato
J.
Caviedes-Voullième
D.
García-Navarro
P.
2016
Rainfall/runoff simulation with 2D full shallow water equations: sensitivity analysis and calibration of infiltration parameters
.
Journal of Hydrology
536
,
496
513
.
FLO-2D Software, I.
2018
Two-Dimensional Flood Routing Model, User Manual
.
Flumen-Group
2010
IBER Modelización bidimensional del flujo en lámina libre en aguas poco profundas, Manual básico de usuario
.
Froehlich
D.
2003
User's Manual for FESWMS FST2DH Two-dimensional Depth-averaged Flow and Sediment Transport Model. volume Report No. FHWA-RD-03-053. Federal Highway Administration, Washington DC, USA
.
García-Navarro
P.
Murillo
J.
Fernández-Pato
J.
Echeverribar
I.
Morales-Hernández
M.
2019
The shallow water equations and their application to realistic cases
.
Environmental Fluid Mechanics
19
,
1235
1252
.
Jaeger
R.
Tondera
K.
Pather
S.
Porter
M.
Jacobs
C.
Tindale
N.
2019b
Flow control in culverts: a performance comparison between inlet and outlet control
.
Water
11
,
1408
.
Liang
D.
Özgen
I.
Hinkelmann
R.
Xiao
Y.
Chen
J. M.
2015
Shallow water simulation of overland flows in idealised catchments
.
Environmental Earth Sciences
74
,
7307
7318
.
Martínez-Aranda
S.
Fernández-Pato
J.
Caviedes-Voullième
D.
García-Palacín
I.
García-Navarro
P.
2018
Towards transient experimental water surfaces: a new benchmark dataset for 2D shallow water solvers
.
Advances in Water Resources
121
,
130
149
.
Murillo
J.
García-Navarro
P.
Burguete
J.
Brufau
R.
2007
The influence of source terms on stability, accuracy and conservation in two-dimensional shallow flow simulation using triangular finite volumes
.
International Journal for Numerical Methods in Fluids
54
,
543
590
.
Nan
Z.
Sheng
J.
Congfang
A.
Weiye
D.
2017
An integrated water-conveyance system based on Web GIS
.
Journal of Hydroinformatics
20
,
668
686
.
https://iwaponline.com/jh/article-pdf/20/3/668/200348/jh0200668.pdf
.
Norman
J.
Houghtalen
R.
Johnston
W.
1985
Hydraulic Design of Highway Culverts
.
Volume Report No. 5 FHWA-IP-85-15 of Hydraulic Design Series
.
Federal Highway Administration
,
Washington, DC
,
USA
.
OpenKinect
2012
Openkinect open-source project
.
Pereira
M. J. M. G.
Fernandes
L. F. S.
Macário
E. M. B.
Gaspar
S. M.
Pinto
J. G.
2015
Climate change impacts in the design of drainage systems: case study of Portugal
.
Journal of Irrigation and Drainage Engineering
141
,
05014009
.
Schall
J.
Thompson
P.
Zerges
S.
Kilgore
R.
Morris
J.
2012
Hydraulic Design of Highway Culverts 3rd ed.. Volume 5 of Hydraulic Design Series
.
U.S. Department of Transportation, Federal Highway Administration
,
Washington DC
,
USA
.
Vreugdenhil
C.
1994
Numerical Methods for Shallow Water Flow
.
Kluwer Academic Publishers
,
Dordrecht
,
The Netherlands
.
Yu
C.
Duan
J.
2017
Simulation of surface runoff using hydrodynamic model
.
Journal of Hydrologic Engineering
22
,
04017006
.