Analysis of the performance of different culvert boundary conditions in 2D shallow flow models

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.

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.

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 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.

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.

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.

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 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.

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 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 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.

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.

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.

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).

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.

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.

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 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.

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.

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.

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.

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.

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.

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.

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.

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.

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