## 2 Materials and Methods  
### 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].  
\[\frac{\partial\rho}{\partial t}+\
abla\cdot(\rho U)=0 \tag{1}\]  
\[\frac{\partial pU}{\partial t}+\
abla\cdot(\rho U\otimes U)=\
abla\cdot(\tau- \rho\cdot\overline{u\otimes u})+S_{M} \tag{2}\]  
where \(U\) is flow velocity of fluid and \(t\) is time step. \(\rho\) is density of fluid and \(S_{M}\) is momentum source. \(\tau\) is the shear stress. Hence, the quantity of \(\
abla\cdot(\rho U\otimes U)\) can be represented by \(\frac{\partial}{\partial x_{i}}(\rho U_{i}U_{j})\).  
In this study, \(k\)-\(\epsilon\) model was applied and the equation is as follows:  
\[\frac{\partial(\rho k)}{\partial t}+\
abla\cdot(\rho Uk)=\
abla\cdot\left[ \left(\mu+\frac{\mu_{t}}{\sigma_{k}}\right)\
abla k\right]+P_{k}-\rho\epsilon \tag{3}\]  
where \(k\) is kinetic energy and \(\epsilon\) is kinetic energy dissipation rate. \(\mu_{T}\) is turbulence viscosity (=\(C_{\mu}\rho k^{2}/\epsilon\)). \(C_{\mu}\) is \(k\)-\(\epsilon\) model constant (=0.09). \(\sigma_{k}\) is the turbulence model constant for the \(k\) equation.  
\[\frac{\partial(\rho\epsilon)}{\partial t}+\
abla\cdot(\rho U\epsilon)=\
abla \cdot\left[\left(\mu+\frac{\mu_{t}}{\sigma_{\epsilon}}\right)\
abla\epsilon \right]+\frac{\epsilon}{k}(C_{\epsilon 1}P_{k}-C_{\epsilon 2}\rho\epsilon) \tag{4}\]where \(C_{e1}\) and \(C_{e2}\) are the Reynolds stress model constants, and \(P_{k}\) is turbulence production due to viscous and buoyancy forces.  
\[P_{k}=\mu_{T}\
abla U\cdot\left(\
abla U+\
abla U^{t}\right)-\frac{2}{3}\
abla \cdot\mathrm{U}(3\mu_{T}\
abla\cdot U+\rho k)+P_{kb} \tag{5}\]  
In the case of the full buoyancy model, the definition of \(P_{kb}\) is as follows:  
\[P_{kb}=-\frac{\mu_{T}}{\rho\sigma_{\alpha}}g\cdot\
abla\rho \tag{6}\]  
where \(g\) is the gravitational constant (=9.81 m/s\({}^{2}\)). \(\sigma_{\alpha}\) is constant.