Large-scale flooding causes widespread disaster, and harmful pollutant concentration in water following flood affects public safety and the environment. In this study, a numerical model for solving the 2D shallow water equations and the solute transport equation is proposed to simulate overland flood and pollutant transport caused by floods. The present model is verified by comparing the predictions with the analytical solutions and simulation results; sufficiently high computational accuracy is achieved. The model is also used to simulate flood inundation and pollution spread in the area of Hun and Taizi Lane (HTL) in China due to river dike breaches; the results show that the coupling model has excellent performance for simulating the flooding process and the temporal and spatial distribution of pollutants in urban or rural areas. We use remote sensing techniques to acquire the land coverage in the area of HTL based on Landsat TM satellites. The impacts of changed land use on mitigation of flooding waves and pollutant spread are investigated; the results indicate that the land cover changes have an obvious influence on the evolution process of flood waves and pollutant transport in the study areas, where the transport of pollutants is very dynamic during flood inundation in HTL area. Furthermore, the motion of pollutants considering anisotropic diffusion is more reasonable than that due to isotropic dispersion in simulating pollutant transport associated with the flood in urban or farmland environments.

1. Introduction

Floods are triggered by the rapid rise of the water level after heavy rainfall; they can cause the river to overflow the embankment, thereby increasing the flood risk in downstream and low-lying areas. Usually, urban environments and rural farmland are vulnerable to flood damage due to the high population density and large areas of crops. Up to now, many areas worldwide have suffered heavy rains and flood disasters. For example, in July 1982, Nagasaki in Japan suffered a massive flood after heavy rain, resulting in property losses of as much as $2.5 billion. In July 1995, the Hun River and Taizi River basins located in Liaoning Province in China experienced large flood events; multiple levees burst in both rivers and 7.6 × 105 hm2 of farmland were submerged by the flood following a dike break; the economic losses reached ¥34.4 billion. The Brisbane River in Australia flooded after continuous rainfall in January 2011, resulting in a severe flood disaster with large impacts on urban residents and the infrastructure [1]. Flooding events in Europe caused losses of about 4.9 billion euros every year from 2000 to 2012 [2].

Numerous studies on flood inundation in a way can minimize the mitigation of negative consequences of flood events by predicting the movement characteristics of floods in [3, 4]. Several different types of models have been developed for flood hazard simulation and prediction in high-risk zones based on the rich floodplain data sources and advanced computing hardware and methods [57]. In earlier studies, a series of simplified hydrodynamic and hydrological models are developed and applied for simulating overland flow processes in the urban [8]. However, it is difficult in these simplified models to obtain accurate estimates of the water depths and velocities of the flood waves in urban areas. In fact, floods tend to have rapid velocity and strong destructive power, so accurate simulations of the water depth and flow velocity of floods in low-lying areas near the river are needed. Therefore, many two-dimensional (2D) models based on the solution to the 2D shallow water equations (SWEs) have become popular due to their computational accuracy and efficiency in flood simulations [1, 9, 10]. Particularly, the Godunov-type SWEs models have exhibited great potential for modeling transient overland flows and flash floods due to their shock-capturing capability [11]. Some researchers carried out a series of discussions of a hydraulic analysis of the impact of land use on the run-off regime in a retention area and consequently on flood wave propagation [12, 13]. The impact of land cover on the propagation of flood waves is important to consider the time and spatial characteristics of the riparian and floodplain areas [14]. Similar studies are carried out to investigate the impact of land uses on flood wave propagation in riparian areas and floodplains; this study also applied constant values of the roughness coefficient [1517]. Vegetation in floodplains significantly reduces water velocities and increases the retention capacity due to levee breaches for the flooding scenarios in the studied area [18].

During the disastrous flood events after intense rainstorms and following dike breaches, in addition to flood mechanical damages, there is also high contamination of floodwater and suspended sediments with inorganic and organic compounds. These pollutants are particularly harmful to public health and the environment and result in deterioration of the water quality [19]. For example, the sewage treatment plants or the treatment ponds that contain pollutants are destroyed by a torrential flood; the released pollutants will rapidly migrate and spread with the floodwater, resulting in significant environmental pollution and damage to densely populated areas. To the best of our current knowledge, an increasing number of numerical studies have been considered to facilitate flood risk management associated with polluted floodwater in urban areas and then to predict the fate of point-source pollutants transported by flood waves in small-scale regions [1, 11, 20, 21]. The mixing and dispersion processes of wastewater surcharging from urban drainage systems within an overland flow are analyzed systematically based on a 2D flood model [22]. A solute transport model based on the Lagrangian particle method is presented to predict the potential contamination paths of pollutants from drainage water in urban areas during a pluvial flood event [23]. Furthermore, a coupled 1D/2D model for hydrodynamic and pollutant transport is developed to carry out the sewer overflow simulation; the results indicate that the model is suitable for simulating the inundation and pollutant transport processes in urban areas [24]. A transport model based on the Lagrangian particle method is proposed to study the potential contamination paths of solutes in drainage water in an urban area during a pluvial flood event [25]. The impact of land uses on nitrogen transport is also investigated to explore how nitrogen transport responds to land-use change and its effects on aquatic habitats in the catchments [26]. A new integrated modeling framework was developed by coupling the Storm Water Management Model (SWMM) to the Cellular-Automata Fast Flood Evaluation model to predict the inundation depth and pollutant transport caused by surcharges over the capacity of the drainage network [27]. In view of the above summary, the land-use change is an important factor for impacting the flood propagation in residential areas, how it interacts, and to what degree it affects flood inundation and pollutant transport that are poorly understood. Under this condition, detailed studies should be conducted to assess the impact of the land-use change on flood inundation and pollutant dispersal in large-scale rural areas.

In this paper, we propose a coupled 2D model to simulate flood disasters and the associated pollutant transport in flood-prone areas. The study is organized as follows: in Section 2, the formulations of the governing equations and a description of the numerical model are presented. In Section 3, the integrated model is validated using several case studies, and the simulated results are compared with the analytical solutions and the measurement data to verify the computation accuracy of the coupling model. Moreover, the study tends to discuss the impact of changed land covers on flood propagation and the pollutant mixing and transport in the areas of Hun and Taizi Lane, China. Lastly, the main conclusions are drawn.

2. Material and Methods

The modeling framework used in this study entails two-dimensional shallow water equations coupled with pollutant mass conservation. The finite volume method with unstructured triangular cells is used to discretize the shallow water equations and pollutant transport equations. The Roe approximate Riemann solution is used to compute the water momentum flux on the interfaces of the triangular mesh. A high-resolution scheme for pollutant transport and diffusion is adopted to deal with the advection term, and a flux limiter is used to reduce artificial diffusion and oscillation. Furthermore, a second-order scheme is adopted to discretize the diffusion term of the solute transport equation.

2.1. Governing Equations

Under the assumption of hydrostatic pressure, the vertical velocity is negligible; a set of SWEs for depth-averaged 2D overland flows over a horizontal plane are deduced [2830]. The SWEs, which are coupled to the depth-averaged solute transport equation, are written in a compact form:where t is the time, U represents the vectors of the conservative variables, F and G are the convection fluxes in the x- and y-directions, respectively, Fd and Gd are the diffusion fluxes in the x- and y-directions, respectively, and S is source term. These parameters are expressed as follows:where h is the water depth, u and are the depth-averaged velocity in the x-direction and y-direction, respectively, is the gravity acceleration, C is the depth-averaged pollutant concentration, and Kxx, Kxy, Kyx, and Kyy are the depth-averaged mixing diffusion coefficients, which can be written aswhere c is Chèzy’s coefficient; and are the dimensionless coefficients of longitudinal dispersion and turbulent diffusion, respectively. In the absence of field measurements, Falconer recommended typical values of  = 13.0 and  = 1.2 [31]. and are the bottom slopes in the x-direction and y-direction, , , and Zb is the bed elevation. SC indicates the concentration added to the source term per control unit, which is obtained from the wastewater disposal through outfalls. D is the inflow discharge. and are the bed friction stress values in the x-direction and y-direction; they are defined aswhere n is Manning’s roughness coefficient at the bottom.

2.2. Numerical Discretization

A Godunov-type finite volume method with an unstructured triangular mesh is used for the numerical discretization and solving of the governing equations [32]. The study domain is first divided into a set of control cells, and then an unstructured triangular mesh is created. An explicit time-marching scheme is used in the calculation model. We use the Green formula for integration with the hydrodynamic equation in equation (1); the discrete value of the next time step is obtained from the previous time step as follows:where and are the average values of the control volumes in the current time step and the last time step, respectively, and Ai is the area of the ith mesh element. The solute transport equation (equation (2)) is also integrated using the Green formula; it is expressed aswhere Qij represents the flux of the boundary j in the control volume I; is the pollutant diffusion flux.

The physical variables are the same in every triangular mesh element due to the average depth integration in the equation; a series of piecewise functions are developed in the solution domain. Hence, the unstructured cells exhibit discontinuity at the boundary of each control unit, resulting in a local Riemann problem. The Riemann problem at the cell interface can be solved by using various Riemann approximations for assessing the interface convection fluxes. In this study, the Roe scheme is used to solve the local Riemann problem [21].

The governing equation of the conservative solute transport is the advection-diffusion equation [33]. When convection is dominant in pollutant transport, the problems of numerical oscillation and high artificial diffusion are common in the first-order upwind scheme. The first-order upwind scheme is improved to increase the accuracy by using the flux limit function [34]. Based on the MUST scheme, an r-factor algorithm using the local element to transfer the upwind information is used to reduce the numerical diffusion error due to the upwind approximation of the advection term in this study [35]. A central difference scheme is adopted to discretize the diffusion term in equation (6), which involves the neighboring points in the center of each element. Consequently, the scheme has second-order accuracy and does not require any restrictions regarding the angles of the triangle meshes.

2.3. Variable Time Step

It is common knowledge that the computational stability for an explicit scheme model is limited by the Courant–Friedrichs–Lewy (CFL) condition and the Peclet (Pe) number [36]; therefore, the variable , which is adaptable to the variability of the solute and hydraulic parameters, is dynamically updated. The time step in the explicit scheme must be satisfied with ; these two dimensionless numbers are defined as follows:where di is the gird’s distance from center to vertices.

3. Simulation Implementation

3.1. Oscillatory Flow in Parabolic Bottom Topography with Friction

In this case study, the oscillation of the free water surface in a parabolic terrain is simulated with bed resistance to verify the ability of the hydrodynamic model for capturing the wet-dry front. The analytical solution of the water level and the bottom topography of the parabolic basin in the oscillatory flow motion are given by [25]. In this experiment, the computational domain is discretized to 18672 triangular grids, the boundary conditions of the study area are fixed wall boundaries, and the initial speed of the study area is 0 m/s. It should be noted that because the boundary of the study area is always in a dry state, the type of boundary condition of the domain does not affect the simulation result. The total simulation time of the water movement is 6000 s, including four and a half cycles. A variable time step is adopted. Figure 1(a) shows the comparison of the numerical and analytical solutions of the water depth at the three gauges (−2750, 0), (−50, 0), and (2750, 0). The measuring points (−2750, 0) and (2750, 0) are located at the junction of the dry and wet areas, and the point (−50, 0) is in a submerged state. The numerical predictions are in perfect agreement with the analytical solutions and no oscillations are detected at the wet-dry interface. The water levels’ numerical results of the model and the analytical solution at times 2000 s and 6000 s are shown in Figure 1(b). A strong agreement is observed between the numerical results and analytical solutions of the free surface profile at different times, and the movement of the intersection between the dry and wet areas is accurately captured, with no signs of spurious distortions. The amplitude of the water surface elevation decreases with an increase in the simulation time. At t = 6000 s, the water level has stabilized at 10 m, the flow is stable, and the wet-dry intersection point is located at x = ±3000 m. The results demonstrate that the proposed method is able to determine the moving boundary of the wet and dry areas.

3.2. Solute Transport Simulation

We validate the performance of the proposed model for solving the advection and diffusion equation of solutes in a constant-speed flow field. A pollutant is released at the source, and the migration of the pollutant is calculated. The test is carried out in a 200 m × 800 m rectangular domain, with the north and south boundaries as wall boundaries, the west boundary as the inflow boundary, and the east boundary as the outflow boundary so that free flow is achieved [1]. The entire domain has a flatbed, and the influence of bottom friction is not considered. The initial conditions are h = 1.0 m, u = 1.0 m/s, and  = 0.0 m/s. The analytical value of the initial pollutant concentration field is given as follows:where C0 = 233.06 kg, and x0 = 0, y0 = 100, and (x0, y0) is where the pollutants are dropped. The mixing coefficients in the x-direction and the y-direction are Kxx = 1.02 m2/s and Kyy = 0.094 m2/s, respectively. In this case, we set the analytical solution concentration at t = 60 s as the initial concentration. In this case study, the computation domain is discretized into 34646 unstructured triangular elements. The simulation time is 600 s, and the time step is variable. Figure 2 depicts the solute concentration profiles over time at y = 100 m. As the pollutant is transported, the concentration decreases from 1 at the initial state to 0.10 at 600 s. The predicted results are in good agreement with the analytical solutions. The temporal and spatial variations of the predicted pollutant concentration are shown in Figure 3. Under the action of the steady velocity in the x-direction, the pollution mass expands rapidly in the flow direction as it is transported downstream. The diffusion coefficient is larger in the x-direction than in the y-direction, and it is observed that the transverse diffusion amplitude is much higher than that in the longitudinal direction. In summary, the results confirm that the solute transport model accurately simulates the diffusion and movement of contaminants in water.

3.3. Simulation Flood Propagation in the HTL

The Liao River basin is one of the largest hydrological systems in China; its watershed area is about 21.9 × 104 km2, and the main river is the Liao River. The Hun River and the Taizi River are two main tributaries of the Liao River, and the two rivers flow southwest in parallel. After the confluence of the Sanchahe in Haicheng City, the river is called the Da Liao River. The Da Liao River empties into the Liao River estuary at Yingkou City and eventually flows into the Bohai Sea. In the section of the Hun River and the Taizi River, the left embankment of the Hun River and the right embankment of the Taizi River form an alley, which is referred to as the HTL (Figure 4). The HTL includes the Xiaobeihe area in Liaoyang County, the Gaotuozi area in Haicheng, the Tangmazhai area, and Wenxiang city. Where the terrain of the HTL is high in the north and low in the south, the maximum difference of terrain elevation is about 11 meters in the whole study domain. The overall terrain is flat and the fluctuations are less dramatic.

Heavy rainfall has a significant impact on the water levels of the Hun River and Taizi River during the flood season. According to historical statistics, flood and waterlogging disasters frequently occurred in the Hun River and Taizi River. A total of 114 floods of varying degrees have occurred in the nearly two hundred years from 1798 to 1997; among them, there are 20 large-scale flood disasters. For example, the city of Majiapuzi near the Taizi River has seen many levee breaches with the heavy flood inundation in 1949 and 1950, and in Daobazi close to the Hun River, a levee break occurred during the great flood of 1995. Large floods have caused significant damage to the urban and rural areas in the HTL. In this study, we use the HTL as a case study and simulate the processes of dike break flood caused by heavy rainfall with a 100 year return period. Daobazi and Majiapuzi are chosen as the dike break sites to simulate a major flood in the HTL; the locations of Daobazi and Majiapuzi are shown in Figure 5. Figure 6 presents the discharge hydrographs of the two dike breach sites. The initial state of the study area is a dry bed, the total duration of the simulation is 120 h, and the initial Manning coefficient is n = 0.025 s/m1/3 in the entire domain. A variable time step is adopted in the simulation to reduce the running time. In shallow water models, an increase in the cell sizes improves the simulation accuracy, but the computational cost increases accordingly. Firstly, we study the flood motion of the dike breach in Majiapuzi station at three different mesh resolutions to determine the most appropriate grid accuracy; five gauges are selected to analyze the processes of flood waves and their location is shown in Figure 5. As shown in Figure 7, the time series of water depth (Figure 7(a)) and the maximum depth (Figure 7(b)) at different gauges corresponding to three mesh resolutions are compared, and the time series of the maximum submerged area (Figure 7(c)) under three resolution grids are also revealed. We can see that the simulated values of water depth and the maximum water depth at the selected gauges almost have no difference under the condition of three grid resolutions, and the submerged area at several gauges is also basically consistent. So we select the mesh accuracy with a total of 9585 unstructured triangular elements, which can complete the calculation at the minimum computational cost and maintain sufficient accuracy. The length of the grid edge around the inflow boundary is small (30∼50 m); the spatial grid size is large in the regions outside of the dike break site to reduce the calculation time. Comparisons of the submerged area are made between the proposed model and the HEC, CUIDES, and FEM models as shown in Figure 8, the results show relatively few differences in the submerged area between these models, and the simulated results from the present model show good computational accuracy. Figure 9(a) shows the inundation range and flow field distribution of the Majiapuzi flood at different times. The floodwater moves from the northeast to the southwest, which is consistent with the terrain of the “Huntai lane” that is higher in the northeast and lower in the southwest. After the dike breach, the submerged area continues to expand, and the floodwater is accumulating downstream. Figure 9(b) shows the flood inundation range and flow field distribution of the Daobazi dike breach at different times. The floodwater moves from the northeast to the southwest. The flood velocity decreases significantly after the inflow discharge has weakened and the flow has stabilized.

4. Results and Discussion

4.1. Flood Simulation Scenarios for Changed Land Uses

In the simulation of overland flow caused by a dike breach of the river, important field parameters which affect the model performance are Manning’s roughness coefficient values in different catchments. HTL is one of the important residential areas and the grain-producing areas in north China, the roles of farmland vegetation and residence area have a crucial impact on the flooding wave processes. It is known to all that remote sensing (RS) data, acquired by optical and microwave sensors mounted on satellites and airplanes, are important alternative sources of information for mapping land use/cover. The objective of the present research is to adopt RS technology to obtain the distribution of residence area and farmland vegetation in the domain of the HTL. The RS images of the HTL area are acquired from Landsat8 Operational Land Imager (OLI), provided by the USGS (https://glovis.usgs.gov/next/). In this computation domain, the emphasis of information extraction is building, waterbody, and farmland vegetation. The normalized difference vegetation index and modified normalized difference water index are used to extract objects. The decision tree classification is used to realize the extraction of buildings and houses, farmland areas, and waterbodies (Figure 10). Land use/cover classes observed in this figure are vegetation land (82.21%), built-up areas (6.71%), and waterbodies (0.92%) in HTL. The value of Manning’s coefficient on the basis of available data on the distribution of land use is set to be 0.13 and 0.05 in the residential district and farmland area, respectively [37, 38]. We select four stations in farmland (F1, F2) and residential areas (R1, R2) to simulate the variation of water depth and flow velocity. As can be seen from Figures 11(a) and 11(b), owing to the rapid flow rate of a flood near the dike break position, the amount of water is temporarily accumulated in a short time at F1 and R1 under the action of farmland and architecture resistance, which increased the water depth at these two stations. The flood arrival time of F2 and R2 gauges is delayed by 2.5 h and 5 h, respectively. Figure 11(c) shows the simulated velocity of flood waves under the condition of different Manning coefficients; the results indicate that the increase of bed resistance attenuates the peak velocity of flood waves and its amplitude decreases by about 40%; the arrival time of the peak velocity is delayed as well. As shown by these values, the presence of the vegetation induced a significant decrease in velocity values. In the later stage of flood evolution, with the increase of floodwater volume, the water depths at the gauges keep constant, where it reaches the state of hydrostatic equilibrium. Generally speaking, the changes of land use in the presence of the residential building, vegetation, and farmland crop not only significantly reduce the effects of the flood wave in the study area but also cause the mitigation of the flood wave and so its effect decreases through the downstream part of the studied area.

We assume that during the evolution of the flood in the Daobazi dike breach, a sewage treatment plant is destroyed by the floodwater in Zhujiajie Village, Wenxiang Town; this allows us to model the transport process of pollutants in a flood of HTL. The solute release point is placed at S1 as shown in Figure 5. The contaminant transports start at the 10th hour after dike break, while the concentration is C = 20 mg/L and the release of pollutants lasts for 10 hours in this case study. In the simulation of the Daobazi dike breach, Figure 12 shows the pollutant transport associated with the overland flow at different times without and with considering the land use types. It is obvious from comparing Figure 12 that the pollutant plume appears a narrower shape using a single Manning’s roughness value in the computation domain compared to that in the situation of considering different land uses, these are all attribute to the variation of overland flow velocity caused by the resistance of farmland and building. After 20 hours at the end of the release, we observe that the highest concentration of the pollutants moves along the movement direction of the flood waves in both cases. The role of advection transport is extremely important compared to diffusion transport during the pollutant process evolution in the flood events; therefore, if the flood waves move fast, the pollutants will quickly spread downstream with the flood waves.

Simulation experiments are also conducted with the pollutant transport subjected to isotropic and anisotropic diffusion driven by the river dike breach flood waves. We set the mixing coefficients εl and εt in equation (3) to 13 and 1.2 for anisotropic diffusion, respectively. For isotropic dispersion, the mixing coefficients are both set to 13. In simulating the pollutant transport following the flood, the effect of farmland and residence resistance is considered; the evolution of the pollutant plume at different times for the anisotropic is compared with the isotropic case in Figure 13. Because of the isotropic dispersion, the pollutant in floodwater spreads faster in a transverse direction than that due to anisotropic diffusion. It is worth mentioning that the locations of the peak pollutant concentration are mainly controlled by the advection. With the increase in transport time, the region of low pollutant concentration will be larger than that under the action of anisotropy diffusion; the concentrations in pollutant centers are significantly low due to the effect of isotropic dispersion. This is caused by the excessive spreading of the pollutant in the transverse direction under isotropic diffusion. The overall development of the pollutant concentration movement and distribution for the anisotropic diffusion is found to be rational.

5. Conclusions

In the present paper, a fully coupled model on solute transport following a flood is proposed to simulate the flow hydrodynamics and the process of pollutant transport in flood-prone districts. Subsequently, the involved model is successfully validated using several tests with satisfactory results, confirming its numerical stability and computational accuracy with shock capture capacity. Moreover, this study shows that the model performs well for a prolonged flood in the HTL area with complex 2D topography during a short time in a large-scale domain. Considering the potential influence of the land-use change on flood events, we focus mainly on the response of flood wave propagation and associated pollutant transport to land cover patterns changes in residential areas. In this simulation, we consider a spatial variation of the roughness using a mesh in which each node has a roughness coefficient depending on the type of land uses; the flood hydraulics are also compared. The results indicate that water depth and flow velocity are sensitive to the variation of the different land-use types; the propagation process of flood waves slows down after considering the resistance of farmland and residential areas, so any considerable change in land use significantly influences the behavior of flood risk during flood events in HTL. In particular, this study also indicates that the pollutant transport is more dynamic following flood inundation in the HTL area, and the land cover plays an important role in the mixing and transport of pollutants in the rapidly varying flood events. In addition, the results of anisotropic diffusion are more reasonable than those of isotropic dispersion in simulating pollutant transport in a real flood study. In general, the present model can quickly predict flood inundation and flood-related water pollution for real-world flood simulations.

In fact, HTL areas have distinct types of plants: grassland, shrubs, forests, or crop. Flood inundation changes from partial submergence to complete submergence for different plants, the resistances caused by plants for different land covers vary with the water depth. One can see that the drag force caused by plants is able to totally describe the flow and resistance in vegetated areas; therefore, this method should be understood and quantified in future work.

Data Availability

All data included in this study are available upon request by contact with the corresponding author.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This work was supported by the National Nature Science Foundation of China (51879028), the National Key R&D Program of China (2019YFC1407704), and the Open Fund of State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology (LP2009).