Three Dimensional Numerical Modeling of Flow and Pollutant Transport in a Flooding Area of 2008 US Midwest Flood

This paper presents the development and application of a three-dimensional numerical model for simulating the flow field and pollutant transport in a flood zone near the confluence of the Mississippi River and Iowa River during the US Midwest Flood in 2008. Due to a prolonged precipitation event, a levee along the Iowa River just upstream of Oakville, Iowa broke, and the small town was completely flooded for a couple of weeks. During this period, the high water level in the flood zone reached about 2.5 meters above the ground, and wind was the major force for the flow circulation. It was observed that some pollutants were leaked from the residential and farming facilities and transported into the flood zone. Leaking of pollutants from these facilities was reported by different news media during the flood and was identified using high resolution satellite imagery. The developed 3D numerical model was first validated using experimental measurements, and then applied to the flood inundated zone in Oakville for simulating the unsteady hydrodynamics and pollutant transport. The simulated pollutant distributions were generally in good agreement with the observed data obtained from satellite imagery.


Introduction
The Midwestern United States is one of four US geographic regions, consisting of 12 states in the United States: Illinois, Indiana, Iowa, Kansas, Michigan, Minnesota, Missouri, Nebraska, North Dakota, Ohio, South Dakota, and Wisconsin.This region has a long history of flooding, with major floods occurring several times over the recent century including 1927, 1961, 1993, 2008, and 2011.Minor flooding is a regular occurrence.
Due to a long term precipitation from December 2007 to June 2008, the Midwest had experienced one of the most wettest winter and spring, causing the soil to become so saturated that additional rainfall could quickly became runoff.In June 2008, much of the Midwestern portions received over 300 mm of rainfall as several storm systems sequentially impacted the region.The vast majority of the precipitation across the region was channeled directly into the lakes, rivers, and streams as runoff, resulting in stream flows reaching historic heights across the Midwest, particularly in many areas of Iowa, south-ern Wisconsin, and northern Illinois [1,2].
Huge flooding events were reported along the Mississippi River and its tributaries.A number of rivers overflowed for several weeks at a time and levees breached in numerous locations.Historical record high stream flows occurred in some of the major regional rivers including the Des Moines, Cedar and Wisconsin Rivers.Reported river crests exceeded 500-year levels in some locations.It was reported that about twenty-two levees were breached forcing evacuations and causing extensive damage by June 20, 2008.Several states, including Illinois, Indiana, Iowa, Michigan, Minnesota, Missouri, and Wisconsin, were affected by the flooding events.Navigation along the Mississippi River was affected with several of the lock and dam structures closed.The flood left two dozen people killed, 148 injured, approximately 35,000 -40,000 people evacuated from homes and the damage region-wide was estimated to be in the tens of billions of dollars [1].
A levee of the Iowa River near Mississippi River failed in the evening of June 14, and the water in Iowa River completely flooded the town of Oakville, and thou-sands of acres of farmland were under water.Since the flooded water was controlled by the levees of Iowa River, Mississippi River, and a hilly area of the Mississippi Valley, an inundation zone was formed for several weeks [3,4].As there was no outlet, after the water rose to the highest level, it kept for several days.It was observed that some pollutants, such as sewage, pig farm waste, diesel fuel and farm chemicals were leaked and transported in the water [5].To study the transportation of the pollutant, the inundation zone was treated as a natural lake.
Due to the lack of field measured data, numerical modeling is one of the most effective approaches to study the pollutant transport processes in the flooding zone.Some well-established models, such as ECOMSED developed by Hydroqual, EFDC and WASP developed by the US Environmental Protection Agency (USEPA), CE-QUAL-ICM developed by the US Army Corp of Engineers (USACE), have been used by many researchers to study the flow, water quality and pollutant transport in natural water bodies [6][7][8][9].Dortch et al. [10] applied the CH3D-WES hydrodynamic model, originally developed by USACE, to simulate the flow field and pollutant transport of Lake Pontchartrain due to the large amount of contaminated floodwater induced by Hurricane Katrina being pumped into the lake.
A three dimensional numerical model, CCHE3D, has been developed by the National Center for Computational Hydroscience and Engineering (NCCHE) to simulate the flow and sediment transport in natural water bodies [11].In this study, a 3D model was developed based on CCHE3D to simulate the flow fields and pollutant transport in the Oakville flood inundation zone.After the flooding wave propagation due to the levee breaching, water level rose to the highest position.The water movement in the inundation zone was mainly driven by the wind and the wind-induced flow was the major force for the pollutant transport.This study presents the validation and application of a 3D numerical model for simulating the wind-driven flow and associated pollutant transport in the flooding area of Oakville during the Midwest Flood.Remote sensing techniques using different satellite imagery acquired during the flood event were used to provide some important information, including topography of study area, locations and dimensions of building and farming facilities, and traces of pollutants on the water surface, etc., for the model's inputs and validation.

Model Description
CCHE3D is a finite-element-based unsteady three-dimensional (3D) free surface turbulence model that can be used to simulate turbulent flows with irregular boundaries and free surfaces.This model is based on the 3D Reynolds-averaged Navier-Stokes equations.By applying the Boussinesq assumption, the turbulent stresses are approximated by the eddy viscosity and the strains of the flow.There are several turbulence closure schemes available including: parabolic eddy viscosity model, mixing length model, k- model and nonlinear k- model.This model has been verified against analytical methods and experimental data representing a range of hydrodynamics phenomena [11].
To study the wind-driven flow in the Oakville inundation zone, a parabolic vertical distribution of eddy viscosity was specified to analyze the wind-induced flow.A mass transport model was developed based on CCHE3D hydrodynamic model.The pollutant transport equation was solved using finite element method which was consistent with the numerical method employed in the CC-HE3D model.

Governing Equations
The governing equations of continuity and momentum of the three-dimensional unsteady hydrodynamic model can be written as follows: where u i (i = 1, 2, 3) are the Reynolds-averaged flow velocities (u, v, w) in Cartesian coordinate system (x, y, z); t is the time;  is the water density; p is the pressure;  is the fluid kinematic viscosity; is the Reynolds stress; and f i are the body force terms.
The turbulence Reynolds stresses in Equation ( 2) are approximated according to the Bousinesq's assumption.These stresses are related to the mean velocity gradients and the eddy viscosity coefficients.The horizontal eddy viscosity coefficient is computed using the Smagorinsky scheme [12]: where x and y are the grid lengths in x and y directions, respectively.The parameter  ranges from 0.01 to 0.5.
For wind-driven flow, the vertical eddy viscosity is mainly determined by the wind stresses.Several formulas were presented to calculate the wind-induced eddy viscosity (Section 2.2).
The free surface elevation ( s ) is computed using the following equation: Copyright © 2013 SciRes.AJCC where u s , v s and w s are the surface velocities in directions;  s is the water surface elevation.
The pollutant transport equation can be expressed by: x, y and z

Wind-Induced Eddy Viscosity
In a natural lake or the ille flood inundation zone, wind stress is the important driving force lation.For simulating the wind driven flow, th tion of vertical eddy viscosity is required.A parabolic eddy viscosity distribution was proposed by Tsanis [13] based on the assumption of a double logarithmic velocity profile.The eddy viscosity was expressed as: in which  is the numerical parameter; z b and z are the characteristic lengths determined at bottom respectively; u is the surface shear velocity; H is the water depth.To use this formula, three parameters, , z b and z s have to be specified.For some real cases with very small water depths, it may cause some problems using this formula to calculate eddy viscosity.Koutitas and O'Connor [14] proposed two formulas to calculate the eddy viscosity based on a one-equation turbulence model.Their formulas were:     ,max 1 5 0.5 1 where h,  is the numerical parameter;  z is the nondimensional elevation, 0.25 ,max 0.105 0.3 0.142 . formula was In this study, a new proposed to calculate ed in a laboratory flume with st the wind-induced eddy viscosity based on experimental measurements conduct eady-state wind driven flow reported by Koutitas and O'Connor [14].The vertical distribution of eddy viscosity was measured and used to fit a calculation formula.
In the proposed formula, the form of eddy viscosity was similar to Koutitas and O'Connor's assumption (Equations ( 7) and ( 8)) and expressed as: Based on measured data, a formula was obtained to calculate the vertical eddy viscosity: Figure 1 compares the vertical distributions of eddy viscosity obtained from experimental measurem formulas provided by Tsanis (Equation ( 6)), Kou (E ents and titas quations ( 7) and ( 8)), and the proposed formula (Equation (11)).The proposed formula can be used to calculate the eddy viscosity over the full range of water depth.At water surface ( 1 z H  ), Equations ( 6) and (8) show the eddy viscosity is zero or very closed to zero.The newly proposed formula shows due to the effect of wind shear, the eddy viscosity at the water surface is about 16% of the maximum eddy viscosity ,max t  .This formula was used to calculate the vertical eddy viscosity in this study.

Boundary Conditions
Wind stress is generally the dominant driving force for flow currents in natural lakes.The the free surface are expressed b wind shear stresses at y where a  y co is the air density; U wind and V wind velocit mponents at 10 m elevation in x tions, respectively.Although the drag coeff d ith are the wind and y direcicient C may vary w wind speed [14,15], for simplicity, many researchers assumed the drag coefficient was a constant on the order of 10 −3 [16,17].In this model, C d was set to 0.001, and this value is applicable for simulating the wind driven flow in several natural lakes [18,19].
On the free surface, the normal gradient of pollutant concentration was set to be equal to zero.In the vicinity of the solid walls and bed, the gradients of flow properties were steep due to wall effects.The normal gradient of concentration on the wall was set to be zero.
There are many dry areas in the flooded zone, which represent the buildings of farmer facilities and islands, and wet and dry technique was applied to those areas.In the model, non-slip and the zero normal gradients boundary conditions were set for simulating the flow velocity and pollutant concentration, respectively.as developed based on the finite ahedral (F on-unifo vection.In this model, a convective interpolation function is used for this ined by solving a linear and

Numerical Solution
The numerical model w element method.Each element is a hexahedral with three levels of nine-node quadrilaterals, and the governing equations are discretized using a 27-node hex igure 2).
Staggered grid is adopted in the model.The grid system in the horizontal plane is a structured conformal mesh generated on the boundary of the computational domain.In vertical direction, either uniform or n rm mesh lines are employed.In order to get more accurate results, the non-uniform mesh lines are placed with finer resolution near the wall, bed and free surface.
The unsteady equations are solved using the time marching scheme.A second-order upwinding scheme is adopted to eliminate oscillations due to ad purpose.This function is obta steady convection-diffusion equation analytically over a one-dimensional local element shown in Figure 2.Although there are several other upwinding schemes, such as the first order upwinding, the second order upwinding and Quick scheme, the convective interpolation function is selected in this model due to its simplicity for the implicit time marching scheme [11].The velocity correction method is applied to solve the dynamic pressure and enforce mass conservation.Provisional velocities are solved first without the pressure term, and the final solution of the velocity is obtained by orrecting the provisional velocities with the pressure solution [11].The system of the algebraic equations is solved using the Strongly Implicit Procedure (SIP) method.
Flow fields, including water elevation, horizontal and vertical velocity components, and eddy viscosity parameters were computed by CCHE3D.After getting the effective source terms S, the pollutant concentration distribution can be simulated by solving pollutant transport Equation ( 5) numerically.ing the proposed formula (Equation ( 11)).Near the bottom, the simulation obtained by using the proposed formula gave better predictions.

Model Verification for the Mass Transport
The proposed numerical model was tested against an analytical solution for predicting the concentrations of a non-conservative substance in a hypothetical one-dimensional river flow with constant depth and velocity.At the middle of the river, there is a point source with constant concentration (Figure 4).Under this condition, the concentration of the substance throughout the river can be expressed as: where U is the velocity; C is the substance concentration; D is the x dispersion coefficient; x is the displacement from point O; and K is the decay rate.An analytical solution given by Thomann and Mueller [20] is: where C 0 is the concentration at O point.In this test case, it was assumed, the water depth is 10 m and D is 30 m 2 /s, the values of K are 0, 1.0/day and 2.0/ Figure 5 shows the concentration distributions obtained along the Iowa River just upstream of Oakville failed, and the water in the river rushed through the broken levee and completely flooded the small town.Figure 6 shows the flood zone and the location of the levee breaching.As t r levee and hill nd the high water level in this zone tented to remain for several weeks.It was reported that the length of the damaged levee was about 1100 m, and it was not completely fixed until July 1 [2,3,21].It was estimated from the high water mark, the high water level was about 2.5 meter above the ground [22,23].It was also observed that some oily pollutants were leaked and transported on the water surface [5].To study the pollutant transport in the inundation zone, this flooding area can be treated as a natural lake.odel, it was applied to simu-he flooding area was governed by rive y area, an inundation area was formed a s As this area is relatively flat, the wind is the major driving force of the water flow.To demonstrate the capailities of the developed m b late the wind-induced flow and pollutant transport.
Figure 7 shows the computational domain and topography obtained from 30 m resolution of USGS Digital Elevation Model (DEM).The inundation zone was an enclosed water body with the north-south axis spanning about 16 km, while the shorter east-west axis was about 8 km.There were many residential houses and farm facilities in this flooding zone, which could affect the flow fields.In the numerical model, they were represented with a number of dry nodes.The computational mesh was generated using the CCHE Mesh Generator.In the horizontal plane, the computational domain was represented by a 199 × 300 irregular quadrilateral mesh.In the vertical direction, it was divided into 8 non-uniform layers and mesh lines conformed to the curved and irregular bed bathymetry.The vertical distribution of the non-uniform grid was generated using the flexible and powerful two-direction stretching function with finer spacing near the surface and bottom [24].

Application of Satellite Imagery for Providing Information for Numerical Model
During the flooding events, it is almost impossible to   8) and detect and map the locations and sizes of residential houses and farm facilities (Figure 7).This imagery along with another Landsat 5 TM imagery acquired on June 16, 2008 were used to map the traces of pollutants derived by flood water and their nature as well.
Both Landsat and ALOS imagery were visually inspected to map the traces of the pollutants.The imagery was displayed in true color and false color composite with 4 (red), 3 (green), 2 (blue) band combination (Figure 9).Principal Component Analysis (PCA) was performed and i udies determine mage derived spectral profiles were st the nature of the pollutants [26].to The spectral profiles derived from the visually detected pollutant traces clearly reveals that the pollutants are some kinds of fluid but not water (Figure 10).From the visual inspection it was observed that the pollutants fluenced by winds.The traces were observed on both Landsat and ALOS imagery maintaining the same pattern  were floating on the water surface maintaining trails in-(Figure 9).The shape and pattern of these trails clearly indicate that the materials composing the traces were no quickly soluble in the water.By the available satellite data it was not possible to determine the true nature of the pollutant materials since there was no measured data for the pollutants.However, since during the flood event a large amount of gasoline leakage was reported [5], the detected pollutant traces in satellite imagery were considered to be of oily pollutant.

Flow Simulation
Water flow in the flood zone was primarily wind-driven.Pollutant stripes on the water surface were observed from satellite image acquired at 16:47 hours, June 21, 2008.Those pollutants were apparently leaked from farming facilities, and floated on the water surface.It was observe ind-induced surface flow velocities.To study the flow d by NCCHE for simulating two dimensional depth-averaged unsteady flow, senundation/propagation, levee inundation due to the levee breaching along the Iowa River near Oakville, and the computed results were set as initial conditions for the 3D model simulation.
The water discharge that entered the Oakville flood zone can be estimated based on the measured discharges at several stations located along Mississippi River and Iowa River, including Lock and dam #17, Keithsburg, and Wapello stations (shown in Figure 6).The following formula was used to calculate the discharge entering the flood zone (Q 0 ): where Q i is the flow discharge at Wapello in Iowa River; Q m is the discharge at Keithsburg in Mississippi River; Q 17 is the discharge at Lock and dam #17, about 5 km upstream of the confluence of Mississippi River and Iowa River.
zone, causing the water level rising gradually from June Figure 11 shows the discharge that entered the flood d that their transports were majorly affected by the w 14 to June 24.In CCHE2D model, the discharge was set as inlet boundary for simulating the flood propagation.
Figure 12 shows the simulated processes of flood propagation/inundation at different times.It can be found that field as well as pollutant transport on the water surface, the 3D model is needed.In this study, the developed 3D model was applied to simulate the flow field induced by wind, and the pollutant concentration distributions were simulated by solving Equation ( 5) numerically.
To run the 3D model, the initial flow fields are needed.CCHE2D model was applied to provide initial flow field conditions, including water surface elevation, velocity, shear stress, eddy viscosity, etc., for the 3D model.The CCHE2D model was develope diment transport, flood i breaching and pollutant transport [27,28].It has been successfully applied to simulate the flood events and flood propagation during the Hurricane Katrina and Midwest flood 2008 [29,30].In this study, CCHE2D model was used to simulate the processes of flood propagation/  pt rising until June 21.The simulated results showed that the averaged water depth was about 2.65 m, close to the observed data (2.5 m).The 2D flood simulation results were used as initial conditions for the 3D model flow simulation.
Due to a long period of time submerged under the water, some oily pollutants were leaked from farmer's facilities and transport in the flooding water.Satellite image (Figure 9) showed that there were some pollutant strips observed on water surface on June 21, 2008.However, one day before, there were no pollutant strips in the same area.It can be estimated that the pollutant leak might occur on June 21.In this study, a period from 6:00 to 18:00 hours, June 21, 2008 was chosen for modeling the flow and pollutant transport in the flooding zone.
Figure 13 shows the ctions near the study area during the simulation period, his study, Equation (11) was us simulated processes of almost all the residential houses and farm facilities were inundated within first 6 hours.After one day, the whole area of the tiny town was under the water, and then the water level ke observed wind velocities and dire and the major wind direction is blowing from northwest, which is consistent with the pollutant moving direction.After obtaining the initial condition provided by the 2D model, the unsteady wind-driven flow fields were simulated using the developed 3D model.Based on some pretests, simulated wind-induced flow fields show similar patterns by using Equations ( 6) and (11) to calculate the vertical eddy viscosity.In t ed to calculate the vertical eddy viscosity in the flood zone.
Figures 14(a) and (b) show the simulated velocity vectors at 4:00 pm on the water surface and near the bottom, respectively.In general, wind-shear forced the flow near the surface to move in the direction of the wind and resulted opposite movement in the lower layers.The water depth and surrounding wall boundaries also affected the circulation of wind driven flow.

Pollutant Distribution
Satellite image showed that the oily pollutants leaked from farming facilities drifted with the flow forming many stripes aligned with the wind and surface flow direction (Figures 9 and 15(a)).After obtaining e simuas applied to simulate ing zone.Due to a nts th lated flow fields, the 3D model w the oily pollutant transport in this flood lack of directly measured data, the amount of polluta leaking was estimated from the satellite imagery and they were simulated as passive materials.The density of oily pollutant is lighter than water; it may float on the water surface, and move mainly following the water surface velocity.
Several locations such as A, B1, B2, C, D, E and F, were selected for simulating the pollutant distribution.In the numerical model, it was assumed that the oily pollutant just moving on the water surface due to the buoyancy force.The leaking rates of the oily pollutant from different locations were set as source terms in Equation (5) and calculated by: in which, S i is the leaking rate of the pollutant at location i; M i is the total amount of pollutant at location i leaking into the water; T i is the time period for pollutant at location i leaking into the water; and V i is the volume of surface water layer receiving the pollutant during the T i time period.The values of M i , V i and T i can be estimated based on the satellite imagery (Section 5).
The simulated pollutan tions at the locations of A, B1, B2, C, D, E and F e compared with the satellite imaginary ( Figures (a) an ).The main trend of pollution stripes was well predicted.Some differences between the model results and llite imagery remain due to the uncertainty of curring time and leaking amount of pollutant.In addition, only one wind station data was available for the simulation and it was applied to the entire flood zone.In reality, wind velocity and direction might change locally over the flood zone.

Effect of Vertical Distribution of Computational Mesh on Surface Flow Fields
It was observed that the oily pollutant transports in the flood zone were mainly affected surface flow velocities.To calculate the wind drive acing near the surface and bottom.
and the total w f surrface velocity increases.The differty obtained using Mesh NonUniform2 ) were very In this study, the flow fields computed using 4 meshes with different vertical distributions were compared, and the ratio between the vertical layer depths ater depth from top to bottom were listed in Table 1.
In the uniform mesh, the depth of surface layer was 12.5% of the total water depth, which was same as the other layers.
In other three non-uniform meshes, the depths of surface layers were 8%, 6% and 4% of the total depth, respectively.Figure 16 shows the simulated surface velocities at the locations of Pt1, Pt2, Pt3 and Pt4 (Figure 7) using the four meshes with different depths of surface layers.It can be observed that with the decrease o face depth, the su ences of the veloci (L1 = 6%H) and NonUniform3 (L1 = 4%H e 1.The ratio between the vertical layer depth (Li) and the total depth (H) of the four meshes.small.Therefore, in the study, the Mesh Non-uniform3 was selected for the simulating the flow and pollutant transport in the flood zone.In this study, remote sensing technology was applied to estimate the length of the broken levee, the pollutant leaking period and the pollutant distribution.Based on the pollutant distribution pattern obtained from satellite imagery, it was observed that most of the pollutant shows a strip following the wind direction.Figure 14(a) shows that the wind-induced surface velocity generally follows the wind direction.So the movement of the pollutant is mainly determined by the wind-induced surface flow velocity.This is one of obvious properties of oily pollutant moving on the water surface.

Useful Information
The satellite imagery was also used to estimate the time and location of the pollutant leaking.However, it is quite hard to estimate the amount of pollutant leaking into the flooding water, and the leaking processes.For implification, the amount of pollutant leaking from each s farmer's facility was estimated based on the polluted area shown on the satellite imagery, and th than 1 , it may be broken [31].ince the aces of pollutants detected in the satellite magery did no dy, the oil po in which, A i is the area of polluted area; d i is the thickness of oily pollutant (0.001 m); and  o is the density of oily pollutant (~0.9 kg/m 3 ).After obtaining the polluted area, the volume of surface water receiving the pollutant (V i ) in Equation ( 18) can be calculated by (20) in which, L 1 is the depth of surface water layer.In the model, it was assumed that the pollutant leaking at several locations occurred simultaneously.Based on the field observation and satellite imagery, the pollutants in the studied flood zone were considered as oily material.It was floating on the water surface and mainly driven by the wind-induced surface velocity.The simulation of wind driven flow was tested using several meshes with di stributions, and the mesh with thinner surface and bottom layers (4% of total depth) was selected for the simulating the flow and pollutant transport in the flood zone.The simulated pollutant distributions were generally in good agreement with observations obtained from satellite imagery.This model provides useful information for understanding the flow circulation and pollutant transport in the flooding zone.It may help for developing and implementing an emergency response plan

Figure 1 .
Figure 1.Comparison of vertical eddy viscosity formulas and experimental data.

Figure 2 .
Figure 2. Coordinate system and computational element (in the figure,  ,  and  are axes of a local system).

Figure 7 .
Figure 7. Computational domain.conductfield measurements to acquire flow field and pollutant distribution data in the inundated areas due to inclement weather conditions.Remote sensing has been widely used as a proven tool to map and monitor flood propagation for many years.During a flood event, nearreal-time satellite imagery serves as a very useful management tool for authorities coping with the disaster[25].Although the observed data by satellite imagery are limited on the water surface, over the years it has been successfully used for water quality studies.In this study,

Figure 12 .
Figure 12.The flood propagation at different times.
the water body was divided into 8 non-uniform layers ertically with finer sp v The computational meshes with different vertical distributions may produce different surface flow velocity.
wat r with 1m thickness.The polluted rea was obta ned y mapping the boundary of the oily ant at location i leaking into the water (M i ) c pollutant leaking (6:00 am) to 16:47, which is 10 hours an6.Summary and ConclusionsA three dimensional model was developed and applied to simulate the wind-driven flow and natural waters.The model was first validated agains measured velocity profile of wind-driven flow obtained experimental open channel.Then it simulate the wind-driven flow and pollutant transport in a flood inundation area in Oakville, IA during the 2008 US Midwest flooding event.The remote sensing technology was used to provide the locations of residential house and farmer's facility, and the length of the broken levee.It was also used to identify the pollutant distributions in the computational domain and extract pollution source locations.

Figure 6. The flood zone and the location of the levee brea- ching.
Initial Bed Elevation (m)

Obtained from Satellite Imagery for Model Simulation
The satellite image acquired at 16:47 hours, June 21, 2008, the time period T i in Equation (18) can be simply assumed from t e