Evaluation of Infiltration Rainwater Drainage (IRD) System with Fully 3-D Numerical Simulation Approach

: Water scarcity can mean scarcity in availability due to physical shortage, or scarcity in access due to the failure of institutions to ensure a continuously regular supply or due to a lack of adequate infrastructure. Water scarcity will be exacerbated as rapidly growing urban areas place heavy pressure on water resources. To solve these problems, various solutions have been applied, but a fundamental solution has not been applied. Recently, a researched and developed infiltration rainwater drainage (IRD) system is being applied with consideration of its applicability. In this study, features of surface runoff and infiltration according to various flow patterns were analyzed using a three-dimensional CFD (Computational Fluid Dynamics) model for calculating water flow in the IRD system. To estimate the optimal setup, a permeability test and scaled model simulation were performed. The runoff characteristics of the IRD system with respect to rainfall intensity and duration were analyzed with dimensionless variables. With the prototype model, the drainage characteristics of the IRD system were analyzed over time using the hydrological curves. From the simulated results, it was found that the IRD system analyzed in this study was appropriate in the field by comparative analysis with the existing system based on peak runoff, internal storage, and lag time. Therefore, by applying the IRD system in the future, it is expected that the IRD has benefits, such as delayed lag time, surface runoff decrease, and an attenuation of the peak runoff.


Introduction
Uncertainty in rainfall patterns due to climate change, urbanization, and industrialization, water security in the urban area has been continuously and globally considered as an important issue. In modern urban pavement, bituminous materials such as asphalt and impervious materials, such as concrete, are used, and a continuous increase in the runoff coefficient was expected. Numerous studies have analyzed the increase of surface runoff and various problems resulting from it with statistical analysis results using water level, rainfall distribution, and population density known to increase the risk [1][2][3]. In particular, changes in urban runoff characteristics cause a decrease in the efficiency of existing rainwater drainage systems, which are known to be a representative cause of various urban water circulation problems due to a decrease in the groundwater level. Global warming due to the climate change can bring the anomalies in the precipitation in the world [4]. In particular, the annual change in rainfall in terms of rainfall intensity in South Korea from 1961 to 2003 was studied by Kim et al., (2005). It was found that the rainfall intensity exceeding 20 mm/h gradually increases and the possibility of developing into a heavy rain further increases [5]. In addition, if the climate change and the increase of rainfall intensity continue, the total rainfall on the Korean Peninsula will continue to increase by the end of the 21st century and it was predicted to increase by about 25% [6]. The impervious land area in Korea has continuously increased over the past four decades, and it is now about 8.0%. Further, in Seoul, the biggest city in Korea, more than 54% of the area was impervious, so rainwater drainage efficiency was very low (Figure 1).

Figure 1.
Increase of Impermeable Area in South Korea [7].
Based on the statistical data for the last 7 years, it was also found that the impermeable area ratio of major cities has increased, except for some island areas (Figures 2 and 3).  Therefore, Park and Shin (2014) conducted a research that the maximum outflow reduction effect of 15% of the total outflow volume and 18% of the peak flow could be expected by controlling the impervious pavement of application site [8]. A rainwater drainage system must have the ability to respond to hydraulic, hydrological, and environmental issues, and problems that are compounded by continuous urban development and climate change. The currently installed conventional rainwater drainage system is a method of draining the surface runoff flowing along the ground surface through grating-type rainwater harvesting on the outside of sidewalks and roads. However, its efficiency was rapidly reduced, due to the illegal dumping of garbage, the inflow of floating debris, and the installation of illegal covers in commercial facilities, which is a structural disadvantage ( Figure 4). In order to overcome the disadvantages of the conventional rainwater drainage system, an Infiltration Rainwater Drainage (IRD) system was developed, and various studies of evaluation and improvement were conducted. The IRD system was one which installed pavers with optimized permeability to utilize the porous volume and underground storage space of pavement materials and reduce surface runoff. In general, this system basically consists of a surface layer, a sub-base layer, and a filter layer [9]. The IRD system is known to have the effect of reducing peak discharge and flood level, and also has ecofriendly functions, such as securing groundwater refilling and mitigating the urban heat island phenomenon [10]. Despite these advantages, the IRD system was not yet widely applied. In particular, it was found that permeable pavement had a physical limitation of the drainage efficiency due to contaminants and pore constriction [11]. Because permeable pavement has a relatively low strength limit compared to pavement materials, such as concrete and asphalt, it was limitedly applied to parking lots, sidewalks, and bicycle paths where the load was light [12]. Therefore, this study was conducted to quantitatively evaluate the drainage performance of the IRD system, which is being actively researched and developed recently to solve the water circulation problem in urban areas caused by climate change. Among the existing numerical studies for the evaluation of the IRD system, the studies that applied the one-dimensional numerical analysis method usually proved the effect of reducing peak discharge due to IRD infiltration through flood level analysis in river basins. On the other hand, existing studies using the three-dimensional, numerical simulations have a limitation of calculating the surface runoff and infiltration separately and then merging the results. Consequently, previous simulations cannot calculate the runoff over time accurately. Therefore, we conducted fully 3-D numerical simulations of the permeability test model, a scaled-down model, and a full-scale prototype model of the IRD system, and evaluated its performance using a hydrologic curve. First, the permeability to determine the storage capacity and retention function were calculated by modeling two cases of permeability experiments for the design of the drainage layer of the IRD system. Second, the applicability of the IRD system with respect to various rainfall intensities and rainfall duration was verified. With the prototype model, the drainage performance of the system per unit catchment area according to the design factors (variables) was quantitatively compared and analyzed. Finally, a simulation was conducted using an integrated model to understand the drainage performance of the IRD system for a unit catchment area. In addition, the surface runoff and infiltration discharge were verified during the rainfall duration. The performance of the IRD system was quantitatively analyzed by comparing the simulation results and runoff characteristics of the existing rainwater drainage system. Im et al. (2007) conducted permeability tests using five types of pavement materials such as concrete pavement, asphalt pavement, soil, grassland, and permeable pavement. It has proven its suitability as a packaging material for rainwater runoff reduction facilities by analyzing that it secures a penetration rate of about 93% or more until it occurs below the rainfall intensity [13]. Im et al. (2009) noted that it was difficult to quantitatively represent the capacity of the infiltration-storage system because it showed a large difference depending on the engineering specifications of the soil and the Antecedent Precipitation Index (API). It was analyzed that the infiltration rate and storage were constant when the rainfall intensity was over 200 mm/h by analyzing the accumulated infiltration and storage ratio. In addition, the permeability of the permeable paver and the conventional paver were compared to prove that the permeability was kept constant [14,15] by analyzing the permeability test results. Köhne et al. (2011) performed simulations by combining the two-dimensional numerical model Hydrus-2D and AMAL-GAM to simulate surface runoff and infiltration [16]. Belheine et al. studied, using numerical tri-axial tests, the micro-mechanical properties of the numerical material, which were calibrated in order to match the macroscopic response of the real material. Based on the local behavior law, numerical simulations were carried out under the same conditions as the physical experiments (porosity, boundary conditions, and loading). They expected the fully 3-D numerical simulation can be helpful for further rainfall drainage management [17]. The IRD system was more efficient than the existing runoff reduction facilities because of its increasing infiltration rate. Such a result was a guaranteed application of IRD as part of runoff-reducing facilities. Therefore, the IRD system was expected to be viable for practical application with subsequent research for improvements [18][19][20][21][22][23][24]. Park et al. analyzed the efficiency and flow pattern of the permeable concrete and block pavement and suggested the ratios of permeable concrete and block pavement with experimental analysis [25].
These previous research results were classified as a completely integrated model. As a method of simulating and merging surface runoff and infiltration simulation results, it has a limitation that it cannot present runoff characteristics over time, so integrated modeling was required.

Governing Equations
For numerical simulation, it was necessary to convert the size, position, and shape of an object into grid materials. In general, the most widely used grid generation techniques are the Finite Element Method (FEM) and the Finite Volume Method (FVM). In this study, surface runoff, infiltration, and subsurface flow were simulated simultaneously by applying the FVM method. The numerical modeling was conducted using the Reynolds Averaged Navier-Stokes model (RANS) to simulate the surface runoff and its turbulent components. The porosity of the permeable pavement varied 20-40%, which was very large compared to the commonly applied packaging materials. It was difficult to design and perform simulations assuming that the flow inside the inner porous material is laminar due to the characteristics of the permeable pavement. Therefore, this study used a directional loss model developed based on Darcy's law to simulate the internal flow of a porous material. To simulate surface runoff and infiltration over time, flow analysis of each flow region was performed simultaneously with two analysis targets, and the runoff characteristics over time were simulated. First, a permeability test was simulated for the design of the drainage layer. The optimal permeability, which was applicable to the IRD system, was calculated. In order to represent the drainage performance of the unit drainage system with the real-scale model, the permeability of the pavement and the runoff with respect to the surface slope were classified with respect to the time, and the peak flow rate and the maximum runoff volume were compared. In general, the turbulence model can be expressed as the sum of the mean and fluctuation components of the turbulence. The unsteady Navier-Stokes equation was modified to obtain the following RANS model from the continuity and momentum equations [26].
where U is flow velocity of fluid and t is time step. is density of fluid and is momentum source. is the shear stress. Hence, the quantity of ∇ • ⨂ can be represented by .
In this study, k-model was applied and the equation is as follows: where k is kinetic energy and is kinetic energy dissipation rate. is turbulence viscos- . is the turbulence model constant for the k equation.
where and are the Reynolds stress model constants, and is turbulence production due to viscous and buoyancy forces.
In the case of the full buoyancy model, the definition of is as follows: where is the gravitational constant (=9.81 m/s 2 ). is constant.

Flows in Permeable Materials
Applying a loss model developed based on Darcy's law, it was possible to calculate efficiently by preventing unnecessary grid generation. The form of Darcy's model applied to this simulation is as follows: where is dynamic viscosity, is permeability, and is the empirical loss coefficient (laminar flow: / and turbulent flow: /2).

Model Setup of Geometry Data
CFD simulations can be divided in three main steps, (a) pre-processing (definition of the region of interest, numerical mesh generation, initial and boundary conditions definitions, mathematical modeling and numerical schemes definitions); (b) simulation solving (numerical solution of the transport equations system); and (c) post-processing (treatment and analysis of results). Procedure of the CFD modeling using Ansys CFX is depicted in Figure 5. We created geometry by importing an existing DesignModeler in ANSYS package. Grid size of three-dimensional mesh was setup with 5 mm.

Error Analysis
The method of calculating the mass error (Merror) of the model is as follows.
where Vpre is volume of precipitation and Vsur is volume of surface runoff.

Permeability Test Results from Previous Study
It was possible to calculate the permeability of the permeable covering material with a permeability test. From the similar test results, Figures 6 and 7, we applied two permeability values (Kperm= 1.0 × 10 −9 m 2 and 1.0 × 10 −10 m 2 ) [28]. The model of the test was combined with water body and soil part ( Figure 4). Boundary conditions of permeability test were summarized in Table 1.

Permeability Test
The validation of the model was performed by analyzing mass conservation. The discharge calculated in the upper part of the model and the infiltration mass in the lower part were compared for the two test conditions as shown in Figure 8 and Table 2.  As a result of the simulations of the test, the peak infiltration was calculated to be 9.5 × 10 −3 and 8.8 × 10 −4 m 3 /s and each peak occurs at 0.067 and 0.78 s, respectively. The maximum bottom displacement was simulated as 6.6 × 10 −5 m 3 /s and 6.7 × 10 −5 m 3 /s at 7.25 and 71.3 s, respectively. The total drainage performance was simulated as 2.3 × 10 −2 and 2.4 × 10 −2 m 3 for the bottom drainage, respectively. In the case of Pe-10, about 3.53 % of increase was observed. However, the total drainage time was 200 s and 1000 s, almost five times increased. With the relationship between the amount of drainage and the total drainage time, the permeability of Pe-10 (=1.0 × 10 −10 ), which had relatively excellent retardation effect and internal storage performance, was determined to be an appropriate permeability to the drainage layer.

Model Setup of Scaled Model Simulation
The scaled model was a combination of a permeable layer and a drainage layer. It was designed and performed to identify surface runoff characteristics of the IRD system with respect to rainfall intensity and duration. There were five rainfall intensity cases, 10~200 mm/h, and three rainfall durations, 100, 200, and 600 s, respectively. The simulation conditions are summarized in Table 3 with 0.1 s of computational time step, ∆ = 0.1 s.  -I1D2  200  Sc-I1D6  600  Sc-I3D1  30   100  5  Sc-I3D2  200  Sc-I3D6  600  Sc-I5D1  50   100  3  Sc-I5D2  200  Sc-I5D6  600  Sc-I10D1  100   100  2  Sc-I10D2  200  Sc-I10D6  600  Sc-I20D1  200   100  1  Sc-I20D2  200  Sc-I20D6  600 As for the specifications of the permeable layer, the permeability coefficient of an optimized value as a result of a previous research, conducted based on the samples collected from a field observation [29], was applied. Three layers were setup for the numerical simulation (Figure 9).  Table 4 shows the total amount and maximum values of surface runoff and infiltration. The cases of the mass error of 0.5% or less were used for the analysis. The rainwater drainage through surface runoff occurs predominantly for all cases. Figures 10-12 show the dimensionless runoff for the same rainfall time by dividing the infiltration and surface runoff by the total runoff on the y-axis, and dividing time by the total rainfall duration on the x-axis. They also show the drainage characteristics in the case of the same duration with dimensionless variables. In the case of the scaled model, 90% of the total runoff is surface runoff, because the volume of porous media is not enough to hold infiltrated water under surface. Figures 10-12 show that the dimensionless surface runoff remains constant over time and the dimensionless infiltration increases. Therefore, it was revealed that there is an effect of reducing surface runoff if sufficient porous media to hold infiltrated water.   When the product of the total rainfall and the duration is the same, that is, when the rainfall depth (I × t, mm) from Table 3 (2.78 and 5.56 mm) is the same, the runoff characteristics according to the rainfall intensity of the IRD system were analyzed as shown in Figures 13 and 14. They represent infiltration and surface runoff with dimensionless variables, respectively. In the case of a rainfall depth of 2.78 mm, the infiltration scale increases by up to 16.8% and an average of 10.8% as the rainfall intensity increase. In case of rainfall depth of 5.56 mm, the surface runoff scale decreases by up to 23.7% and 43.4%. This reversal of the infiltration scale was analyzed to be caused by insufficient volume of porous media compared to rainfall. Therefore, installing the IRD system, it is necessary to design appropriately considering the amount of rainfall and the duration of rainfall.

Model Setup of Prototype Simulation
Analysis of the prototype model was performed to quantitatively represent the drainage performance of the IRD system according to the design catchment area. In this model, block drainage discharge, drainage layer discharge, drainage hole, surface runoff, and infiltration were analyzed with respect to time. For the rainfall intensity applied in this simulation, 210 mm/h, the 200-year frequency of rainfall intensity in the central region of the Korean Peninsula, was applied. The design catchment area of the model was the same as the specifications of the IRD system actually built in Incheon city, South Korea. Figure 15 shows the schematic design of catchment area. For the validation of the model, it was designed based on the domestic constructing regulations [29]. According to MOLIT, in the case of a cross slope 2% is the construction standard, and a maximum of 5% is allowed. In order to simulate the drainage performance according to the transverse slope and permeability of the drainage layer and pavement, simulations were performed by applying permeable blocks of different permeability and compared to the results of application of existing forms. We selected the part from the catchment area in the field as the application model ( Figure 16). With four values of permeability and three cases of slope values including no slope (0%, 3%, and 5%), conditions of numerical simulation were setup as summarized in Table 5. We built the four kinds of drainage system model for the numerical simulation including normal type ( Figure 16).
To improve the accuracy of the simulation of the prototype model, the mass conservation error of each case for the time interval and grid size was analyzed (Table 6). Based on the result that the error was less than 8%, the time interval and grid size were determined as shown in Table 6.    Table 7 shows detailed results of runoff characteristics of the IRD system according to the various structures. The time of occurrence of the maximum infiltration was simulated to be 176 s and 112 s later, respectively. The peak of surface runoff exists at 587 and 600 s. When Pr-Pe4S10 designed as a multiple-structure was applied, the retention effect was demonstrated by slowing the time of maximum surface runoff through relatively fast maximum penetration compared to the cases of control group. It was found that the maximum infiltration flow rate was reduced by 81.21% compared to Pr-Pe0S10, an impervious drainage system. It can be concluded due to these changes in runoff characteristics, that when the IRD system with a multiple structure was applied to cities and watersheds, surface runoff decreases with retention effects by infiltration ( Figure 17).  The hydrograph of the various structure cases showed that the runoff characteristics of the initial rainfall duration of 600 s were significantly changed according to the surface slope. Table 8 shows the runoff characteristics for each case according to the cross-slope. The maximum penetration occurred at 2.1 × 10 −3 m 3 /s, 112 s in the case Pr-Pe4S10 with 0% cross-slope. In addition, the maximum surface runoff reduction rate according to the slope was 69~81%, and the deviation according to the cross slope was about 12%. The cross slope was analyzed as a major design factor of the IRD system. The maximum surface runoff was simulated as 3.93 × 10 −4 m 3 /s, 7.29 × 10 −4 m 3 /s, and 7.99 × 10 −4 m 3 /s, respectively, as the cross slope increased, indicating a proportional relationship between the slope and the maximum surface runoff. In case of total infiltration, 0.923 m 3 , 0.742 m 3 , and 0.756 m 3 were simulated. Thus, the design of an excessive cross slope may cause an increase in the burden on the drainage system during high-intensity heavy rain ( Table 8). The surface runoff and infiltration of the IRD system with respect to the permeability of the pavement were respectively shown in Figure 15. In the case of the IRD system, it was found that the surface runoff amount was reduced significantly (Figure 18).    Table 9 shows the maximum surface runoff, infiltration, and the timing of each occurrence. The maximum permeability occurs in Pr-Pe2S10, which had the largest permeability, but the decrease in the maximum surface runoff occurred in Pr-Pe4S10, which had the least permeability. In the case of Pr-Pe3Sl0, the time at which the maximum surface runoff occurs increases by 20 s, so it can be concluded that the urban water circulation problem can be solved if the delay time was designed with the permeable design. Figure  16 shows the surface runoff and infiltration of permeable rainwater drainage systems with respect to the permeability of the pavement. When the IRD system was applied, it was simulated that the surface runoff was reduced significantly (Figure 19).  Based on the maximum surface runoff and rainfall duration, the runoff of each model was shown in dimensionless form ( Figure 20). As a result of the dimensionless comparison, it was possible to prevent the performance degradation of drainage facilities due to rapid surface runoff by applying the IRD system. In addition, it was found that the penetration pattern varies greatly depending on the design factor. The infiltration rate and surface runoff rate of each drainage system for total rainfall were calculated as shown in Table 10. Up to 71% of the precipitation applied to the design catchment area can be drained through infiltration with the IRD system.  The detailed runoff, internal storage, and retention rates of each drainage system are compared as summarized in Table 11. In the case of Pr-PE1S10 designed with a simple drainage layer structure, the retention rate was simulated as 94.23%. The single drainage layer design was a relatively unfavorable condition for rainwater drainage compared to the complex structure design. As a result of simulation according to the permeability of Time/Rainfall duration the pavement when designing multiple structures, the highest total amount of infiltration was calculated as 1.11 m 3 for Pr-Pe2Sl0 due to the high permeability of the pavement. As a result of simulation of Pr-Pe4Sl0, Pr-Pe4Sl3, and Pr-Pe4Sl5, the infiltration amount decreased with the increase of the transverse slope. As shown in Table 9, this was identified as a phenomenon caused by an increase in the surface runoff with an increase of the crossslope. As a result of the overall simulation, the displacement was plotted with respect to time. As a result of the simulation, it was found that the outflow characteristics changes according to the design factors of the IRD system. The maximum inflow rate reduction effect of more than 50% was found when the IRD system was applied ( Figure 21).
The maximum peak flow reduction rate was compared to the existing drainage system. The deviation of peak flow according to the design factor was 13% according to the structure, 10% according to the permeability, and 12% according to the slope. In other words, the deviation due to the structure of the drainage system was the most significant. The IRD system with a single structure drainage layer stores 99.4% of the total rainfall inside due to the permeability of the drainage layer. The reduction rate of peak discharge was shown in Figure 22. Compared to the case of Pr-PE0Sl0 (impervious pavement case), all the cases of the IRD reduce peak flow. The desirable reduction of the peak can be achieved combining cross section slope and permeability. Pr-Pe1Sl0 Pr-Pe2Sl0 Pr-Pe3Sl0 Pr-Pe4Sl0 Pr-Pe4Sl3 Pr-Pe4Sl5 Figure 22. Comparison of IRD system prototype model.

Discussion
From the result of the simulation of the previous permeability test, the drainage layer with a permeability of 10 −10 m 2 was chosen as the optimal permeability value. The scaled model was built for the IRD system test and its applicability was verified by applying various rainfall intensities and rainfall durations. As the rainfall duration increased, the surface runoff remained relatively constant with an average of 0.941, 0.942, and 0.943. Under the same rainfall depth condition, the infiltration was reversed when the volume of porous media was not sufficient. In the results of the scaled model simulation, the dimensionless values of surface runoff were almost constant at about 0.94 despite the change in rainfall duration. In the case where the product of total rainfall and duration is the same, that is, where the rainfall depth is the same, the runoff characteristics according to the rainfall intensity of the IRD system were identified. That is, in the case of a rainfall depth of 27.6 mm, the penetration result increased by 16.8% on average and 10.8% on average as rainfall intensity increased. For a rainfall depth of 55.5 mm, the magnitude of surface runoff decreased by up to 23.7% and up to 43.4%. As such, it is judged that the result that the infiltration result was calculated in the opposite way was because the infiltration area was not sufficiently secured compared to the amount of rainfall. Therefore, when constructing the IRD system, it is necessary to design the application site considering the rainfall amount and rainfall duration.
Prototype modeling cases (Pr-Pe2Sl0, Pr-Pe3Sl0, and Pr-Pe4Sl0) were simulated with respect to the permeability of the permeable block, respectively. As a result of the simulation, the total storage rates were simulated to be 88.28, 31.37, and 94.78%, respectively. The permeability of the permeable block of Pr-Pe3Sl0 was the second biggest one among three of permeability variation cases. Although Pr-Pe3Sl0 had the second biggest permeability, it had the least retention rate among the three cases. Through the analysis of the relationship of permeability-total drainage time, in the case of Pr-Pe3S10, rainwater stagnated inside the permeable block and in the drainage layer was effectively drained after rainfall. This design can maximize the retention effect of the infiltration type rainwater drainage system that collects and drains rainwater slowly by using the internal pore internal storage function during rainfall. If the IRD system was designed with the test Peak outflow reduction ratio condition of Pr-Pe3Sl0 (permeable block with no cross slope), infiltration was equivalent to 62% of the total rainfall and 69% was continuously drained up to about 1,200 s under a heavy rain condition. It was concluded that the IRD system can be reduce the burden on the municipal drainage and the quantitative evaluation with 3-D integrated numerical modeling method suggested in this study.

Conclusions
Previous numerical studies have analyzed behavior of surface runoff and infiltration in the permeable pavements, but few studies used only 3-D numerical simulation or coupled modeling after separately analysis. This study was conducted to quantitatively identify the drainage performance of the IRD system with 3-D numerical simulation. Variables such as peak flow rate, peak flow rate reduction rate, internal storage amount, and lag time were analyzed. Compared to the existing system, the relative advantage and applicability were evaluated with a 3-D CFD modeling.
The performance of the IRD system could be evaluated with the procedure suggested in this study. The 3-D CFD modeling, which includes surface runoff and internal infiltration, is expected to be applicable to the optimum design with respect to the dimension and the characteristics of the drainage area and rainfall pattern. The construction materials, such as blocks and sand, used for the drainage system, will vary with respect to construction conditions. While physical modeling will provide accurate evaluation for the various conditions, it takes time and effort to perform the physical modeling. The CFD modeling calibrated with physical model will further improve the accuracy and efficiency of the evaluation the IRD system under the various conditions. With HPC resources, this approach and procedure are very useful for the evaluation of rainfall drainage in urban area.