Porous Shallow Water Modeling for Urban Floods in the Zhoushan City, China

Typhoon-induced intense rainfall and urban flooding have endangered the city of Zhoushan every year, urging efficient and accurate flooding prediction. Here, two models (the classical shallow water model that approximates complex buildings by locally refined meshes, and the porous shallow water model that adopts the concept of porosity) are developed and compared for the city of Zhoushan. Specifically, in the porous shallow water model, the building effects on flow storage and conveyance are modeled by the volumetric and edge porosities for each grid, and those on flow resistance are considered by adding extra drag in the flow momentum. Both models are developed under the framework of finite volume method using unstructured triangular grids, along with the Harten-Lax-van Leer-Contact (HLLC) approximate Riemann solver for flux computation and a flexible dry-wet treatment that guarantee model accuracy in dealing with complex flow regimes and topography. The pluvial flooding is simulated during the Super Typhoon Lekima in a 46 km2 mountain-bounded urban area, where efficient and accurate flooding prediction is challenged by local complex building geometry and mountainous topography. It is shown that the computed water depth and flow velocity of the two models agree with each other quite well. For a 2.8-day prediction, the computational cost is 120 min for the porous model using 12 cores of the Intel(R) Xeon(R) Platinum 8173M CPU @ 2.00 GHz processor, whereas it is as high as 17,154 min for the classical shallow water model. It indicates a speed-up of 143 times and sufficient pre-warning time by using the porous shallow water model, without appreciable loss in the quantitative accuracy.


INTRODUCTION
Worldwide, urbanization develops rapidly over the past decades. The number of urban inhabitants increases from 33.35% of the world population to 55.27% in 2018(United Nations Population Division, 2018, and is expected to 70% in 2050 (Gross, 2016). In China, the urbanization rate increases from 36.22% in 2000 to 54.77% in 2014 (Zhang et al., 2016). Especially for the most economic areas along the east coast, the urbanization rate is about 80% or more. It is suggested that urbanization has substantial impacts on the hydrological processes, including rain-island effects (Min et al., 2011), reduction of land-surface permeability and natural water storage capacity, increase of surface runoff and vulnerability of drainage systems etc. Combined with the climate change that the frequency and intensity of extreme storms and rainfall events increase considerably in space and time with global warming (IPCC, 2014), urban flooding becomes one of the most devastating hazards for the human society. It causes comprehensive damage to people, infrastructure, environment and economy (Huang et al., 2017;Kreibich et al., 2019). For example, in China, more than 100 cities suffer from urban flooding every year since 2006 and over 100 million citizens are evolved (Zhang et al., 2016). Urban flooding may be initiated by heavy rainfall events (pluvial flooding), unexpected dam-break flow, levee overtopping along coasts or major rivers. In practice, it is unlikely to prevent the occurrence of urban floods, but it is possible to reduce flood risks and implement mitigation measures. In this regards, urban flood modeling is a paramount tool to depict the flood prone areas and inundation processes, which are necessary to support emergency management and identify most effective mitigation measures.
In this study, we focus on two-dimensional (2D) urban flood modeling. There have been a number of recent 2D modeling studies on this issue Hunter et al., 2008;Sanders et al., 2008;Liu et al., 2015;Bellos and Tsakiris, 2016;Leandro et al., 2016;Löwe et al., 2017;Xing et al., 2019;Cheng et al., 2020;Liu et al., 2020). Liu et al. (2015) developed a 2D cellular automata model for simulating rainfall-induced urban floods. Bellos and Tsakiris (2016) used a 2D fully shallow water model with two types of infiltration modes for flood simulation in small-scale catchments. Leandro et al. (2016) and Löwe et al. (2017) combined the dual-drainage modeling and hydro-inundation model for urban floods. Liu et al. (2020) coupled the hydrological and hydrodynamic module into an integrated urban flood model and analyze the influences of underlying surfaces. However, the complex hydrological and hydrodynamic processes and fine-resolution demands around buildings of arbitrary geometry often cause computational overhead and thus large-scale simulation (e.g., regional, catchment or city scales) by fully 2D shallow water equations is still a major challenge (Vacondio et al., 2014;Sanders and Schubert, 2019;Xing et al., 2019;Liu et al., 2020). Though simplified diffusive models have been used to reduce the computational burden in conjunction with parallel computations or GPU programming (Hunter et al., 2008;Chen et al., 2012;Leandro et al., 2014;Leandro et al., 2016;Costabile et al., 2017), it is found that the diffusive model might show poor prediction around buildings (Costabile et al., 2017) and could be even less computationally effective than the dynamic model for high resolution meshes (Hunter et al., 2008;Neal et al., 2012). Moreover, simple grid coarsening of structured meshes degrades the building resolution and thus leads to large computational errors Brown et al., 2007).
In order to increase computational efficiency while maintaining adequate accuracy, porosity sub-grid models have been devised and applied for urban flood modeling (Defina, 2000;Guinot and Soares-Frazão, 2006;Sanders et al., 2008;Soares-Frazão et al., 2008;Schubert and Sanders, 2012;Kim et al., 2015;Özgen et al., 2016;Guinot et al., 2017;Sanders and Schubert, 2019;Varra et al., 2020). Such porositybased models enable using relatively coarse computational cells while preserving to some extent the detailed topographic information at a subgrid-scale by means of isotropic or anisotropic porosity parameters. Their applications to experimental-scale modeling show substantial improvement in saving computational cost (e.g., an order of 10 2 ∼10 3 times accelerating) while obtaining the same order of accuracy in water depth simulation (Sanders et al., 2008;Soares-Frazão et al., 2008;Cea and Vázquez-Cendón, 2010;Kim et al., 2015;Özgen et al., 2016;Guinot et al., 2017). Recently, porous shallow water models are attempted to use for field-scale urban flood modeling and show promising results (Kim et al., 2014;Guinot et al., 2017;Sanders and Schubert, 2019). Kim et al. (2014) simulated the Malpasset dam-break flood (France) by a porous shallow water model based on different mesh types. Guinot et al. (2017) used the dual integral porosity model to simulate a hypothetical levee breach flood into the neighborhood of West Sacramento (CA). Using Message Passing Interface (MPI) directives, Sanders and Schubert (2019) applied the porous shallow water model (Sanders et al., 2008) to several largescale cases such as the Baldwin Hills dam-break and Los Angeles region extreme flooding etc.
In this paper, two models have been developed and compared for the city of Zhoushan (China) which is endangered by Typhoon-induced urban flooding every year. One is the classical shallow water model that approximates complex buildings by locally refined grids, the other is the porous shallow water model that adopts the concept of porosity for parameterizing building effects. Both models are developed under the framework of finite volume method using unstructured triangular grids, along with the HLLC approximate Riemann solver for flux computation and a flexible dry-wet treatment that guarantee model accuracy in dealing with complex flow regimes and topography. The pluvial flooding is simulated during the Super Typhoon Lekima in a 46 km 2 mountain-bounded urban area, where efficient and accurate flooding prediction is challenged by complex building geometry and mountainous topography. In addition, the impacts of infiltration processes related to different vegetation types are also discussed.

Governing Equations
The governing equations of the 2D porous shallow water model are composed of mass and momentum conservation equations in x and y horizontal directions (Sanders et al., 2008). The integral form is y ) hu(un x + vn y ) + gh 2 n x /2 hv(un x + vn y ) + gh 2 n y /2 where Ω and zΩ control volume and boundary of an arbitrary cell, respectively ( Figure 1A); zΩ p interface between buildings Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 687311 and voids; r and r p distance measured counter-clockwise around zΩ and zΩ p , respectively; t time; U vector of conservative variables; E n outward normal vector of flux variables; n [n x , n y ] represents the unit outward normal vector of cell edge; s vector of source terms; sp vector of interfacial term representing fluid pressure along the interface zΩ p ; h water depth; u, v depth-averaged velocities in x and y directions; g acceleration of gravity; q rainfall intensity per unit area; f infiltration rate per unit area; s fx and s fy friction slopes in x and y directions estimated by Eqs (3a,b); s bx and s by bed slopes in x and y directions defined by Eqs (4a,b).
where n is the Manning roughness coefficient and z b bed elevation. The building effects on storage and conveyance in the porous shallow water model are represented by the binary density function i′ i′ x, y 1 η′ x, y > z b x, y 0 else where η′ is water level. Note that when i′ is set to 1 for all the cells and neglecting the interfacial term, Eq. 1 becomes the classical shallow water equation. Accordingly, the volumetric porosity ϕ that corresponds to the fraction of cell volume occupied by voids and the areal porosity ψ that corresponds to the fraction of the plane occupied by voids are calculated by (Özgen et al., 2016) where z 0 is the minimum bed elevation in the cell.
The interfacial term s p can be split into a stationary and a nonstationary components by (Sanders et al., 2008) zΩ p where the first and second terms on the right hand side of Eq. 7 are the stationary and non-stationary components, respectively; m [m x , m y ] is the unit normal vector along zΩ p with an inward direction toward the building and m′ [0, m x , m y ]; h| η 0 is the hypothetical still water depth; C b D is the drag coefficient related to the building effects; u [u, v] is the velocity vector and u′ [0, u, v]; V u 2 + v 2 √ . By substituting Eqs (6a,b) and Eq. 7 into Eq. 1, the integral form of the porous shallow water equations can be rewritten as

Numerical Methods
Based on the numerical methods in Hu et al. (2019), both of the porous and classical shallow water equations are solved under the framework of finite volume method using unstructured triangular grids. For simplicity, only the details for the porous shallow water model are elaborated.

Time Integration and Spatial Discretization
The time integration of Eq. 8 over Δt for cell i reads Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 687311 where U n i and U n+1 i represent the vectors of conservative variables at time step n and n + 1, respectively; Δt is time step restrained by the Courant-Friendrichs-Lewy (CFL) condition with courant number of 0.45; The subscript i and j refer to cell i and the j-th edge of cell i, respectively; A i is the area of cell i; Δr j is the length of the j-th edge of cell i. Other variables and parameters are the same as the above unless otherwise specified.

Computation of Flux, Source Terms and Dry-Wet Fronts
For the computation of advection fluxes E n at the interface of adjacent cells, the high-resolution Harten-Latex-van Leer-Contact (HLLC) approximate Riemann solver which is able to automatically capture shockwaves and discontinuities (Toro, 2001) is adopted. For the estimation of bed slope and friction, the slope flux method by Hou et al. (2013) and the splitting pointimplicit method by Liang and Marche (2009) are used, respectively. The dry and wet cells are automatically judged and computed by the threshold water depth of 1 × 10 −6 m. Details are referred to Hu et al. (2019).

Computation of Interfacial Term s p Stationary Component
For the stationary component (i.e., the last term on the right hand side (RHS) of Eq. 9), it is assumed that the hydrostatic fluid pressure along zΩ p is balanced by the bed slope and the hydrostatic fluid pressure along zΩ (Sanders et al., 2008). So we may write (10) where n′ [0, n x , n y ] and s′ b [0, s bx , s by ]. The hydrostatic fluid pressure along zΩ in Eq. 10 can be evaluated by where h j is the water depth at the mid-point of the j-th edge of cell i, which is estimated by Eq. 12 in order to ensure non-negative value (Audusse et al., 2004).
where η 0i is the static water level approximated by the water level at the cell center; Z bj is the maximum bed elevation of the two cell centers located on the two sides of the interface j. The discretization of the bed slope term in Eq. 10 follows the slope flux method by Hou et al. (2013), where h i and Z bi is the water depth and bed elevation at the cell center; Z bj min(Z bj , η 0i ).

Non-Stationary Component
For the non-stationary component, the approach used for vegetation resistance modeling (Nepf, 1999) is extended to estimate the building drag force. Considering the anisotropy property, the non-starionary building drag force per unit area can be written as where K is the number of building edges in a cell; the subscript i is the cell index, the subscript k is the building edge index; C 0 Dk is the drag force coefficient of the k-th building edge; E uk [0, l k m x u, l k m y v]; l k is the length of the k-th building edge. In case the building edge is in the sheltered back face of the building, C 0 Dk is set to zero. Other variables are the same as specified above.

Estimation of Infiltration Effects
The infiltration processes of the land surface are described by the Horton′s method (Akan, 1992), which is where f p,t is the infiltration capacity at time t by the Horton′s method; f 0 and f ∞ are the infiltration capacity at t 0 and for the saturated status, respectively; k d is the decaying rate of the infiltration capacity; f (t) is the actual infiltration rate at time is the total infiltration amount till time t following the Horton′s method; F(t) t 0 f (t)dt is the actual total infiltration amount till time t. Other parameters are the same as specified above.

STUDY AREAS
The Zhoushan archipelago, which is composed of thousand islands, is located at the conjunction of the Yangtze River Estuary and Hangzhou Bay in the East China Sea ( Figure 2A). The Zhoushan city is in the main island of the archipelago, which has an area of 468.7 km 2 and is the fourth biggest island in China. Under the influences of subtropical marine monsoon climate, storms and typhoons occur frequently in summer and autumn and cause heavy rainfall in this area. The rainfall is mostly concentrated in the south-western region (Dinghai district) where the city center is located.
In the period of 1960-2009, there are 187 typhoon events that have affected Zhoushan and 52 events (27.8%) induce heavy rainfall (e.g., daily rainfall >50 mm) (Wang et al., 2011). Based on the historical data, there are 42 and 15 events characterized by the total rainfall higher than 100 and 200 mm, respectively. In the past 2 decades, the typhoon-induced rainstorm is intensified in both frequency and magnitude due to the global warming. Moreover, the lack of big rivers to accommodate strong rainfall flow and the dramatic increase of impermeable surface due to rapid urbanization aggravate the urban flood risks. Consequently, large areas of farmlands and urban districts are inundated frequently during typhoon-induced rainstorms. In addition, the extreme rainfall may also cause flash floods, landslides and debris flow in the mountainous areas. These result Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 687311 in catastrophic destroy in infrastructures, ecology, environment, social economy and even death.
In the present simulation, we focus on the rainfall-induced surface runoff in a catchment of Dinghai district (Figure 2A). The catchment is 46.28 km 2 with a perimeter of about 31.3 km, where dense urban areas and mountains are included.

Rainfall Processes
We consider the rainfall processes when the super Typhoon Lekima passed the Zhoushan City from August 9 th to 11 th in 2019. The variation of rainfall intensity per hour at the Chen′ao  Reservoir (Figure 2A) is shown in Figure 3, for which the data is available online at http://www.zjsq.net.cn:8010/webwarn/. The rainfall mainly occurs in the first 30 h featured by a strong unsteady and undulated pattern. The maximum rainfall intensity is 41 mm/h at 14 h. As no more rainfall data is openaccess for this catchment, we assume the rainfall is evenly distributed in space and the rainfall process at the Chen′ao Reservoir is used for all computational cells.

Building Treatment
In the catchment, different building treatments are employed for the dense urban areas and the isolated houses, respectively. The former is composed of many residential units bounded by streets or borders (e.g., relatively large red polygons in the lower part of Figure 2A). Due to the lack of high resolution data for the details inside the residential units, we consider each residential unit as an integral by neglecting the detailed information of buildings. In comparison, the latter is characterized by a much smaller size and scattered distributions (e.g., small red dots in Figure 2A). The impacts of an individual house on the hydrodynamic processes are spatial and temporal limited, so we make a simplification for saving the computational costs. For an isolated house with an area larger than 3,000 m 2 , we consider its effects by describing its border through polygons. Otherwise, the effects of isolated houses are neglected.

Infiltration Effects of Different Vegetation Types
The infiltration effects of the underlying surface are determined by many factors, such as vegetation types, soil properties, and rainfall intensity etc. (Cen et al., 1997;Liu et al., 2014;Fernández-Pato et al., 2016;Xing et al., 2019). Here, we focus on the vegetation-related infiltration effects. There is a diversity of vegetation in the catchment including grass, crops, vegetables and trees like poplar, bayberry and grapefruit, etc. In the modeling, we categorize them into three types by grass, bush and arbor, respectively. Their spatial distributions are illustrated in Figure 2B.
According to Tao et al. (2017), the influences of different vegetation on the infiltration processes are manifested by the infiltration capacity (f 0 , f ∞ ) and the decaying rate k d ( Table 1). For each computational cell, the parameters of the comprehensive infiltration capacity and the decaying rate are the weighted average of each vegetation type based on the distribution density of each vegetation in the cell.

Numerical Cases, Grids and Parameters
Three cases are designed to investigate the pros and cons of the porosity model and the infiltration effects on the modeling results. Case 1 is the reference case using the porous shallow water modeling (PSWM) without infiltration effects. Its comparison with Case 2 using classical shallow water modeling (CSWM) shows the magnitude of accuracy and efficiency of the porosity method. Case 3 is similar to Case 1 but with infiltration effects.
As shown in Figures 1B,C, the CSWM directly excludes the buildings from the computational domain and refines the grids around the buildings. This is the most delicate way to model the surface water flow along the streets. Yet, it is also time consuming for large-scale simulation. To save the computational cost while maintaining adequate accuracy, the PSWM uses relatively coarse grids to simulate the surface water flow while the building effects of storage and conveyance are considered by the volumetric and areal porosity parameters. Thus, the water flowing area of a coarse grid in the porosity model is equivalent to a random polygon bounded by the edges of the grid and the buildings inside. In this case, for the grids which are entirely inside the building areas, no computation is done; only for the grids which cover the streets (some grids may include small parts of the buildings), the computation is conducted. The total water flowing areas in the two models are therefore the same though the total areas of the effective computational domains in the two models have some differences due to the imperfect overlap of the edges of buildings and coarse grids. In the simulation, the total grid number for PSWM (Cases 1 and 3) is 22,397 with the grid size ranging from 43 to 82 m; it is 161,665 for CSWM (Case 2) with the grid size between 12 and 30 m.
Along the borders of the catchment, the free-outflow boundary condition is implemented. As no big rivers in the catchment, we neglect the river channels for simplicity while the reservoirs in the mountains are reserved. The pluvial flood is initiated by the rainfall of the Typhoon Lekima on an initially dry bed and the simulation lasts for 66 h. Actually, the Manning roughness coefficient n in estimating the bed friction effects (Eq. 3a, b) may be varied for different land surfaces, e.g., the road made of concrete may have different roughness from the vegetated areas. However, for the real urban flood modeling, the Manning roughness coefficient is usually set as constant to avoid uncertainties raising from the empirical value setting for different land surface materials (Schubert and Sanders, 2012;Kim et al., 2014;Sanders and Schubert, 2019). Thus, we use a constant Manning roughness coefficient n 0.02 for both PSWM and CSWM. The empirical parameter C 0 Dk for estimating the drag force is usually considered as a constant (Sanders et al., 2008;Schubert and Sanders, 2012;Kim et al., 2014;Guinot et al., 2017). Its value ranges between 0 and 2 for different experiments and field-scale simulations, and there have been no consensus on which value is more appropriate. In the present modeling, the effects of different buildings on the drag force are accounted for by the vector E uk in Eq. 14 while the empirical parameter C 0 Dk is assumed to be two based on the suggestions by Munson et al. (2006) for square forms in 2D flows. Extra tests (not shown) also indicate that the computed water depth and velocity in PSWM are not so sensitive with different values of C 0 Dk (e.g., 0, 1, 2), but they show better agreement with CSWM in dense urban areas when C 0 Dk 2 is used. All the simulations are carried out using 12

Comparison Between PSWM and CSWM
For the simulation of 66 h pluvial flood in the concerned catchment, the wall clock time of the CSWM and PSWM is 17154 and 120 min, respectively. It indicates that the porosity method may speed up the simulation by about 143 times saving 99% of the computational time. The computations of the two models generally agree with each other quite well, but discrepancies in specific locations are still visible. The spatial and temporal properties of water depth and flow velocity and their comparisons between the two models are shown as following.

Water Depth
The spatial distributions of water depth are qualitatively the same in the two models (Figures 4, 5). Among the most parts of the mountains, the water depth is the lowest. The surface shallow flow runs downward along hillslopes resulting in deeper water in the valleys. The three reservoirs in the mountains retaining the rainfall and surface runoff show a large water depth. At the interface of the mountains and the urbanized areas, the water depth shows a remarkable increase possibly due to the blocking effects of the residential units and the sudden decrease of bed slope. This is exemplified by the large water depth just outside the urbanized areas at the foot of the mountains (e.g., yellow and red parts on the left of Figures 4, 5). In addition, the lowland areas without significant urbanization also show an obvious retaining effect on the rainfall and surface flows (e.g., light blue and green parts on the right of Figures 4, 5). Among the residential units, the water depth along the streets is mostly between those in the valleys and lowland areas. The differences of water depth among those streets depend on the width and directions of the streets, arrangements of the residential units, and changes of land-surface elevation, etc. For the time evolution of the surface flow, the water depth shows an increasing trend with the enhancement of rainfall intensity (Figures 4, 5). Note that there exists a time lag between the rainfall peak (t 14 h) and the maximum water depth (t 15 h) in the catchment (see Figures 3-5). This is mainly because that it takes time for the surface runoff to accumulate and the change of water depth could be delayed. Furthermore, Figure 6 illustrates the magnitudes of water depth increase from 12 to 15 h for the two models showing similar trends. The largest depth increase during the rainfall processes is in the streets of the dense urban areas especially for those near the main path connecting the interface (e.g., red color in Figure 6). In addition, the interface area that is featured by deep water (e.g., orange color on the left of Figure 6) and the partly-blocked lowlying areas (e.g., yellow-green color on the right of Figure 6) also show a remarkable increase. In the mountainous areas the depth increase is rather small (e.g., dark blue in Figure 6) except for the reservoirs.
To get a further insight into the details of the two models, some discrepancies are still existed (Figures 4, 5). On the mountainous areas, the calculated water depth in the PSWM is not as smooth as that in the CSWM. For the former, small puddles are scattered in this area mainly due to the low resolution of the complex topography (e.g., locally steep slopes and delicate undulations) when using large grids. In the dense urban areas, the flow width and depth along the streets are a bit larger in the PSWM. This is probably related to the imperfect overlap of the grids and borders of the residential units  Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 687311 8 when using large grids and the approximation errors at extremely low porosity locations ( Figure 1B).
As the observational data are not easily available, we select 12 points (Figure 2A) to examine the model performance in different region types. Points 1-4 lie in the dense urban areas to represent the street flow, Points 5-6 at low-lying areas partly surrounded by isolated houses and residential units, Points 7-9 at the interface between mountains and urbanized regions, and Points 10-12 in the mountainous areas. The computed water depth variations show distinct characteristics among different regional types (Figure 7). For the highland areas (e.g., P10 and P11), the surface flow is directly determined by the rainfall intensity, so the dry bed becomes wet immediately after it starts raining. The surface water is shallow with remarkable depth fluctuations in response to the sharp changes of rainfall intensity. For other low-lying areas (e.g., other points), the water depth is affected by both rainfall intensity and the convergence of surface flow, so there is a delay for the depth increase and the maximum depth is much larger. Though sharp increase of water depth is also observed, the temporal variation is much smoother compared to that in the highlands. For the comparison of the two models, good agreement is achieved at most points though visible discrepancies are still existed in streets between residential units (e.g., P1-P3) and steepsloped hills (e.g., P11) due to the complex layouts and geometries of buildings and the low grid resolution for steep slopes.

Velocity
Though the velocity computation of PSWM has a lower resolution compared to that of CSWM, the general trends of the two models are very similar (Compared Figure 8 with Figure 9). In the mountainous areas, low velocity occurs on the hillslopes where the surface flow is very shallow. Due to the accumulation of surface flow in the valleys, the flow energy is strengthened and the velocity becomes higher. This results in that the maximum flow velocity is distributed along those valleys. Along the interface of the mountains and the urbanized areas, the velocity depends on the blocking extent and flow mobility. Under strong blocking effects, large amount of water is retained locally and the flow is decelerated dramatically (e.g., the low velocity area with deepest water on the left of Figures 8, 9). In case there is other space available for the flow routing, for example, moving along the interface, the velocity might be as high as that in the valleys. Inside the urbanized areas, the velocity along streets could be similar or even higher than that at the interface. Its spatial distribution depends on the width and directions of the streets, alignments of residential units and bed elevations, etc. In general, the streets featuring deeper water show higher velocity, which bears similarity to valleys. Like the evolution of water depth, the flow velocity in the catchment is also increased with rainfall intensity but characterized by a delayed response.
At the specific points, the flow velocity is more sensitive to the rainfall intensity than water depth ( Figure 10). It is shown by the strong fluctuations of velocity with relatively large magnitude at heavy rainfalls and a smooth trend with small magnitude at light rainfalls (e.g., after 30 h). The velocity could be very small either in very shallow waters (e.g., P1, P4, P10) or in the retained deep waters (e.g., P5, P8, P12). Besides, the magnitude of velocity is also largely affected by local topography. For instance, a steep slope may introduce large flow velocity in shallow waters (e.g., P11). In Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 687311 general, the computed velocity of the two models ( Figure 10) is in line with each other. However, because the velocity computation has a higher requirement on the grid resolution, the velocity differences of the two models are a bit larger than the water depth.
To further compare the direction and magnitude of velocity computation in the two models, the spatial distributions of velocity vectors at t 14 h are shown in Figure 11 (red refers to PSWM, black refers to CSWM). Along the streets in the residential units ( Figure 11A) and around the isolated houses ( Figure 11B), the vectors are generally in good agreements. This implies that the PSWM could adequately capture the main features of flow movement in such areas dominated by the interactions of buildings and streets. In addition, the vectors in the mountains also agree well between the two models ( Figure 11C) showing acceptable accuracy of PSWM with large grids in the steep-sloped mountainous areas.

Impacts of Land Surface Infiltration
The impacts of land surface infiltration are investigated for PSWM by the comparison of Cases 1 (no infiltration) and 3 (with infiltration). The temporal evolution and spatial distribution of depth and velocity are almost unchanged when the infiltration effects are considered, whereas both magnitudes reduce markedly as shown by the computed results at t 15 h in Figure 12. Figure 13 shows the reduction in magnitude and percentage for depth, respectively. On the aspect of magnitude, the largest reduction is observed in the areas of relatively large water depth, shown by dark blue (i.e., the interface between mountains and urbanized areas, and the wide streets connecting the main path to the interface) on the left of Figure 13A. In terms of percentage, the most significant reduction occurs mainly in the streets of the dense urban areas ( Figure 13B). These findings are inconsistent with the spatial distribution of vegetation types. According to the spatial distribution, the vegetation on the mountains is dominated by bush and anchor ( Figure 2B), which have relatively large infiltration capacities, and thus the reduction directly related to infiltration should be larger than that in the urbanized areas where grass is dominant. However, the largest reduction is in the grass-dominated areas. This implies that the infiltration processes not only have a local effect based on the vegetation type, but also extend its spatial influences through surface flow networking. The surface flow in the urbanized areas may have two sources: one is the rainfall, the other is the converging water flow from the mountainous areas. Yet the relatively strong infiltration on the mountains causes less surface flow to go downward to the mountain-city interface and the streets inside, so the convergence contribution is reduced. In case the role of converging flow reduction exceeds that of the infiltration differences between different vegetation types (e.g., between grass and bush/arbor), the water depth decrease could be considerable and this might be the reason for the largest reduction in the urbanized areas.  Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 687311 11 Figure 14 shows the comparison of temporal variations of water depth at the specific points for Cases 1 and 3. By considering the infiltration effects in Case 3, a delay is observed for the start of water depth increase (i.e., land surface from dry to wet) at all points when compared to Case 1 without infiltration. This delay is the longest in the mountainous areas (e.g., P10 and P11) and the shortest in the streets of urbanized areas (e.g., P1-P4). Other places are in-between (e.g., P5-P9). Regarding the reduction magnitude,  Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 687311 the infiltration is mainly effective during the high rainfall intensity periods before t 30 h. After that the difference of water depth with and without infiltration becomes minor. For the flow velocity, the time lag is also existed at all points for Case 3 with infiltration. However, except for the obvious velocity decrease in the streets near the main path areas connecting the interface, other places do not show a universe infiltration effect on flow velocity.

CONCLUSION
Two models (i.e., the classical shallow water model and the porous shallow water model) have been developed and compared for urban flood modeling in the city of Zhoushan (China). The pluvial flooding induced by the super Typhoon Lekima in 2019 in a 46 km 2 mountain-bounded urban area is simulated. The main conclusions include: The computation of water depth and flow velocity of the two models agree with each other quite well. It is shown that the largest water depth is prone to occur in the interface between mountains and dense urban areas where surface flow is dramatically decelerated by building blocking and abrupt flattening in terrains. Moreover, the streets in urban areas connecting the main paths to the interface and the partly-blocked low-lying regions may also encounter large water depth. The highest flow velocity occurs in the steep-sloped valleys among mountains due to the surface flow converging from surrounding hillslopes.
The infiltration processes not only have a local effect based on the vegetation type, but also extend its spatial influences through surface flow networking. Thus, the magnitude of water depth  Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 687311 reduction may not be consistent with the infiltration capacity of vegetation if the networking effect is stronger than the infiltration difference related to local vegetation type. For a 2.8 days prediction, the computational cost is 120 min for the porous shallow water model using 12 cores of the Intel(R) Xeon(R) Platinum 8173 M CPU @ 2.00 GHz processor, whereas it is as high as 17,154 min for the classical shallow water model. It indicates a speed-up of 143 times and sufficient pre-warning time by using the porous shallow water model, without appreciable loss in the quantitative accuracy.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
WL contributed to the idea initiation, result analysis, article writing and supervised BL to develop the porous shallow water model. BL developed the porous shallow water model and did the simulations. PH contributed to the idea initiation and result analysis. ZH contributed to the result analysis. JZ prepared the rainfall, topography and vegetation data. All authors contributed to the article and approved the submitted version.