Advances in Water Resources Numerical simulation of ﬂoods from multiple sources using an adaptive anisotropic unstructured mesh method

The coincidence of two or more extreme events (precipitation and storm surge, for example) may lead to severe ﬂoods in coastal cities. It is important to develop powerful numerical tools for improved ﬂooding predictions (especially over a wide range of spatial scales - metres to many kilometres) and assessment of joint inﬂuence of extreme events. Various numerical models have been developed to perform high-resolution ﬂood simulations in urban areas. However, the use of high-resolution meshes across the whole computational domain may lead to a high computational burden. More recently, an adaptive isotropic unstructured mesh technique has been ﬁrst introduced to urban ﬂooding simulations and applied to a simple ﬂooding event observed as a result of ﬂow exceeding the capacity of the culvert during the period of prolonged or heavy rainfall. Over existing adaptive mesh reﬁnement methods (AMR, locally nested static mesh methods), this adaptive unstructured mesh technique can dynamically modify (both, coarsening and reﬁning the mesh) and adapt the mesh to achieve a desired precision, thus better capturing transient and complex ﬂow dynamics as the ﬂow evolves. In this work, the above adaptive mesh ﬂooding model based on 2D shallow water equations (named as Floodity) has been further developed by introducing (1) an anisotropic dynamic mesh optimization technique (anisotropic-DMO); (2) multiple ﬂooding sources (extreme rainfall and sea-level events); and (3) a unique combination of anisotropic-DMO and high-resolution Digital Terrain Model (DTM) data. It has been applied to a densely urbanized area within Greve, Denmark. Results from MIKE 21 FM are utilized to validate our model. To assess uncertainties in model predictions, sensitivity of ﬂooding results to extreme sea levels, rainfall and mesh resolution has been undertaken. The use of anisotropic-DMO enables us to capture high resolution topographic features (buildings, rivers and streets) only where and when is needed, thus providing improved accurate ﬂooding prediction while reducing the computational cost. It also allows us to better capture the evolving ﬂow features (wetting-drying fronts).


a b s t r a c t
The coincidence of two or more extreme events (precipitation and storm surge, for example) may lead to severe floods in coastal cities. It is important to develop powerful numerical tools for improved flooding predictions (especially over a wide range of spatial scales -metres to many kilometres) and assessment of joint influence of extreme events. Various numerical models have been developed to perform high-resolution flood simulations in urban areas. However, the use of high-resolution meshes across the whole computational domain may lead to a high computational burden. More recently, an adaptive isotropic unstructured mesh technique has been first introduced to urban flooding simulations and applied to a simple flooding event observed as a result of flow exceeding the capacity of the culvert during the period of prolonged or heavy rainfall. Over existing adaptive mesh refinement methods (AMR, locally nested static mesh methods), this adaptive unstructured mesh technique can dynamically modify (both, coarsening and refining the mesh) and adapt the mesh to achieve a desired precision, thus better capturing transient and complex flow dynamics as the flow evolves.
In this work, the above adaptive mesh flooding model based on 2D shallow water equations (named as Floodity) has been further developed by introducing (1) an anisotropic dynamic mesh optimization technique (anisotropic-DMO); (2) multiple flooding sources (extreme rainfall and sea-level events); and (3) a unique combination of anisotropic-DMO and high-resolution Digital Terrain Model (DTM) data. It has been applied to a densely urbanized area within Greve, Denmark. Results from MIKE 21 FM are utilized to validate our model. To assess uncertainties in model predictions, sensitivity of flooding results to extreme sea levels, rainfall and mesh resolution has been undertaken. The use of anisotropic-DMO enables us to capture high resolution topographic features (buildings, rivers and streets) only where and when is needed, thus providing improved accurate flooding prediction while reducing the computational cost. It also allows us to better capture the evolving flow features (wetting-drying fronts).

Introduction
Over the last two decades, the risk of urban flooding in heavily populated coastal regions has increased and is expected to increase further, mainly due to the urbanization and climate change. Urban flooding in coastal regions could be caused by a single source (heavy rainfall, high sea levels or storms), or several sources acting in combination ( Dawson et al., 2008 ). Due to this increasing high risk, the combined effect and the joint probability of multiple extreme floods are gaining importance in flooding simulations ( Lian et al., 2013 ). Improving the predictive capabilities in such cases is critical for populated areas, especially cities. It is therefore important to develop an efficient and accurate numerical model for studying floods caused by several concurrent hazards in coastal cities. bility ( Sto. Domingo et al., 2010 ). Recently various efforts to compute overland flows by solving the shallow-water equations have been made. These studies have simulated overland flows under extreme and unsteady rainfall conditions and spatially constant infiltration rates have been taken into account ( Bellos and Tsakiris, 2016;Nguyen et al., 2016;Singh et al., 2014 ). However, due to the complexity and uncertainty of flood modelling, efficient simulation of flooding at highresolution terrain remains a significant challenge in hydrologic and hydraulic studies. For efficient and accurate flood inundation modelling, numerous methods have been developed, such as grid coarsening methods ( Hartnack et al., 2009 ), cellular automata approach ( Dottori and Todini, 2011 ), and speeding-up strategies such as parallel processing ( Hu and Song, 2018;Sanders et al., 2010;Zhang et al., 2014 ). Nguyen et al. (2016) have proposed a coupled model called HiResFlood-UCI, which combines the hydrological model HL-RDHM and the hydrodynamic BreZo model while ensuring a bare minimum computational cost. HiResFlood-UCI uses HL-RDHM as a rainfall-runoff generator and BreZo as the hydrological routing scheme. This model has been successfully applied to a catchment of the Illinois river in USA.
More recently, a new 2D double control-volume finite element (DCV-FEM) flood model has been developed together with the adaptive unstructured mesh technology and validated by a simple flooding case . In comparison to adaptive mesh refinement (AMR) (a fine structured mesh nested within a coarse mesh) technique ( Berger and Colella, 1989 ), the DMO technique is able to adapt the mesh optimally in time and space in response to the evolving flow features, thus providing sufficient mesh resolution where and when it is required. In this work, we have further developed this adaptive unstructured mesh shallow water model with anisotropic considerations for modelling urban floods from multiple sources (rainfall and storm surge). The implicitscheme has been adopted for solving the shallow water equations and applications to urban floods caused by multiple sources (rains and sea levels). In the DCV-FEM scheme, the velocity components are discretised FE-wise, while the pressure/free surface height is discretised CV-wise .
In this paper, a DCV-FEM adaptive mesh urban flooding model based on 2D shallow water equations named as Floodity , has been further developed for modelling the concurrent flooding and has been successfully applied to a 2.2 km × 1.7 km densely urbanized area within Greve, Denmark. This is the first time to apply the anisotropic-DMO method to simulate urban floods caused by multiple sources based on high resolution Digital Terrain Model (DTM) data. Model validation has been performed in comparison with results from MIKE 21 FM. Sensitivity of the extreme sea levels, rainfall and mesh resolution to the flood volume has been explored to assess uncertainties in model predictions.
The structure of this paper is as follows: The governing equations and anisotropic-DMO method are introduced in Section 2 . Section 3 demonstrates the application of Floodity to simulate flood events within Greve, Denmark. It describes the study site and geospatial data used for modelling, and specifies boundary conditions and parameter settings. Section 4 presents and discusses the Floodity results by comparing with results from MIKE 21 FM. Finally, some conclusions are given in Section 5 .

Methodology
In this work, we have adopted the element pair P 1 DG-P 1 CV  ) (a modification of linear discontinuous velocity and continuous pressure representations) for 2D shallow water simulations. In DCV-FEM scheme used here, the pressure (free surface height) is discretized CV-wise rather than FE-wise in the classic CV-FEM ( Forsyth et al., 1989;Geiger et al., 2004;Gomes et al., 2017;Jackson et al., 2015;Matthai et al., 2007 ). The DCV-FEM provides significant improvements in the quality of the pressure matrix that can be solved efficiently even on highly anisotropic elements. We also use flux limiting of the free surface height based on the Normalized Variable Diagram (NVD) approach of Leonard (1988) .
The Discontinuous Galerkin (DG) method is used for the discretization of the momentum equation. The DG approach is very powerful as it has a natural dissipation associated with it. To robustly stabilize the shallow water momentum equation, a non-linear Petrov-Galerkin scheme ( Farrell and Maddison, 2011 ) is used. It is a mathematically consistent residual scheme and converges to the governing equations as the mesh and time step size are refined. Due to a diffusion term proportional introduced to the residual of the momentum equation, it provides the robustness needed when there are sharp changes in velocity (usually occurs near wetting and drying fronts). A non-linear method is used for the time discretization in which the value of (between 0.5 and 1) is adjusted in space and time. is calculated at each CV face based on the satisfaction of a Total Variational Diminishing (TVD) criterion ( Pavlidis et al., 2016 ). The nonlinear iteration is based on the fixed-point iteration method described in Salinas et al. (2017) . This is important as wetting and drying (due to the non-linear drag and inertia) is highly non-linear.

Momentum and continuity equations
For depth averaged velocity u in non-conservative form, the momentum equation or shallow water is written as: and the continuity equation is written as: where h is the water depth, C f is the volumetric drag coefficient, is the dynamic viscosity, p is the depth averaged pressure, s b is the source term of velocity (unit: m s −2 ), s h is the source term of mass (from rainfall for example, unit: m s −1 ), and p is calculated: and thus where g is the gravitational acceleration, and d is bathymetry (the height deviation from a horizontal and flat plain). The pressure and free surface height are defined CV-wise, contrary to the classical CV-FEM, where they are FE-wise.

Drag coefficient
A commonly used bottom stress parameterization is the Manning-Strickler formulation ( Keulegan, 1938;Rouse, 1961 ): in which n is the unit normal to the bottom surface Γ bottom normal, is the kinematic viscosity, and n m is the Manning coefficient. The formulation for the volumetric drag coefficient C f is: Wetting and drying: Here we use the thin film wetting and drying algorithms, which specify a minimum threshold depth that defines the categories of wet or dry in the model. Importantly, we have introduced the flux limiter to ensure the positive water depth and avoid the physical oscillation. ℎ min = 0 . 01 mm was finally chosen as considering a layer of water of 0.01 mm in the dry areas is physically plausible, while also reducing the non-linear behaviour which may be introduced by a smaller value (e.g. 0.001 mm). In order to enable the CVs to wet and dry freely, we perform an element wise average of the node-wise depths of water within each element. We then use this average in Eq. (6) .
Representation of buildings: The most common methods for simulating the water flow among structures are: (1) blocking-out of the solid area; (2) local elevation rise of the solid area; (3) local increase of roughness of the solid area either via the Manning coefficient increase ( Bellos, 2012 ). The method used here is similar to the local elevation rise. The elevation and shape of buildings are embedded in the realistic topography data. With the use of anisotropic-DMO, the details of buildings can be captured as the water floods around them.

Boundary conditions for the joint flooding events
The boundary conditions need to be set for the water depth h , pressure p or velocity u . We can define either pressure or velocity boundary conditions but not both ( Gresho and Sani, 1998 ). Here we specify the pressure p boundary conditions and then the velocity is calculated from the pressure gradient.
h boundary conditions: For the joint flooding events involved with pluvial and coastal flooding, the water depth h boundary conditions along the coastline include: (1) time series of extreme water level heights e estimated from storm-surge models or the historical data observed; (2) bathymetry (terrain elevation) d along the coastline. Thus, the water depth h is calculated from ℎ = − . p boundary conditions: Once the h boundary conditions are set, the pressure p boundary conditions are specified through Eq. (3) , and the velocity is calculated from the pressure gradient ( Gresho and Sani, 1998 ) and evolve to point into or out of the domain (depending on the dynamics). The boundaries (except the coastline) over the domain can be set up as closed or open based on actual situation.
Rainfall as source term: The rainfall needs to be considered as the source term s h in Eq. (2) . Generally, given a large domain (a large catchment), where the characteristics of the subcatchments are found to be significantly different from each other, the effect of spatial distribution of rainfall is considered. This mean that different subcatchments may receive different amounts of rainfall, namely different rainfall intensity time series are induced according to the region. However, since the computation domain in this study is relatively small, we assume that the rainfall is uniform across the domain. Due to the fact that the effect of rainfall is relatively small compared with that of incoming waves, infiltration is not considered here. An average depth of runoff is then obtained, which contributes to the calculation of the water depth h .

Anisotropic dynamic mesh optimization technique (anisotropic-DMO)
The anisotropic-DMO has the advantage of capturing details of surface and local flows (wetting-drying front) during the process of flooding modelling. The use of anisotropic-DMO can efficiently provide a high mesh resolution where and when it is needed ( Hiester et al., 2011 ). That is, finer meshes are placed only in specific regions where the variations of flow variable solutions are relatively large (e.g. flow around buildings and along the flooding paths), while coarser meshes are used in areas far from these regions, where inundation has not yet occurred for example.
Here, the anisotropic-DMO technique relies on the derivation of appropriate error measures, which dictate how the mesh is to be modified. The required error measure is defined in the form of a metric tensor. The metric is derived from a solution field variable and an error norm based on the interpolation error ( Pain et al., 2001 ) defined to make sure that a desired element length can be obtained while having a required inter-polation error. Thus, the metric ̂ is calculated from where ̂is a scalar constant, and ̂= 1 is used here, ̂is a required interpolation error, and H is the Hessian matrix for a specified field ( Ω) (here, the field of water depth): The desired edge length, h i , in the direction of the i th eigenvector e i , of the metric ̂ , is defined as ℎ = 1∕ √ , i.e., anisotropic as well as inhomogeneous resolution results from a mesh that respects the metric ̂ , where i is the eigenvalue associated with e i ( Piggott et al., 2009;Power et al., 2006 ).
It is advantageous to modify the metric to impose the maximum and minimum element sizes on the mesh, especially for problems that have known high-aspect ratio dynamics or domains. To impose these maximum and minimum constraints directionally, ̂ is modified and defined as where V is a rotation matrix that includes the eigenvectors of the metric ̂ calculated from Eq. (7) and general directions for the maximum and minimum edge length can be introduced through the use of V .
To bound the aspect ratio of elements in physical space, there is a need to limit the ratio of eigenvalues. For example, to limit the aspect ratio of elements to be below a , the eigenvalues are modified as follows ( Pain et al., 2001 ): where  showing the corresponding mesh (bottom row) based on bathymetric data without buildings in scenarios of an individual 2100 upper extreme sea-level event (left column) and a joint event with 100-yr return period rainfall (right column).
in which a is the maximum aspect ratio. h min and h max are the minimum and maximum element sizes. A uniform isotropic element can be transformed to an adapted anisotropic element under the transformation in the physical domain, achieving the desired interpolation error everywhere.  Fig. 1 (b)) simulated using MIKE 21 FM and Floodity with a mesh resolution of 20 m, 10 m, and 5 m respectively, based on bathymetric data without buildings in a scenario of 2-yr return period rainfall and 2100 upper extreme sea-level event.  Fig. 1 (b)) simulated using MIKE 21 FM and Floodity with a mesh resolution of 10 m based on bathymetric data without buildings in a scenario of 2-yr return period rainfall and 2100 upper extreme sea-level event.
The Galerkin interpolation technique ( Farrell and Maddison, 2011 ) is used for interpolating the solutions from the previous mesh onto the newly adapted mesh. The overhead and extra computational cost introduced during the adaptive mesh procedure is very low, around 10% of the total CPU time .

Descriptions of study site and data
To assess the performance of anisotropic-DMO in flooding modelling, the new flooding model has been applied to an urban area located within a 2.2 km × 1.7 km densely urbanized area within Greve, Denmark.  The study area is located in the northeastern part of Greve, Denmark, which covers part of the coastal area (see Fig. 1 (a)). Historical extraordinary flood events which were caused by a series of rain events have occurred in Greve ( Kommune, 2007 ). In addition to extreme rainfall events, Greve is also vulnerable to flood induced by extreme sea-level events along its coast. For example, the most extreme historical flood occurred on 13th October 1760 with a maximum water level of 3.7 m, was caused by a very serious storm surge ( Madsen, 2008 ). Recurrence of any of these flood events would cause numerous damages (e.g. loss of life, direct damages to roads, railways and buildings, indirect damages including loss of income, clean-up cost, turnover loss, cost of illnesses, etc). As a consequence, it is important to develop flooding models to improve the accuracy of flood predictions. In the Greve case study, the digital elevation data provided was quality assured and buildings were incorporated into the DTM (Digital Terrain Model) data with a resolution of 1.6 m ( Fig. 1 (b)), which was detailed enough to describe topographic features (buildings, rivers and streets). Data of different extreme rainfall events as well as extreme sea-level events have been used as the boundary conditions and sources.

Extreme sea-level and rainfall events
Due to the impact of climate change in the next 100 years, future climate change conditions should be taken into account to estimate the future extreme sea-level and rainfall events ( Berbel Roman, 2014 ).
To forecast storm surges, a hydrodynamic model has been calibrated against historical events during 2000-2015. The capability of this built hydrodynamic model to forecast storm surges was evaluated against historical storm event during 2010-2017. An example of a real time forecast issued during the storm "Gorm " on 29th November 2015 can be seen in Fig. 2 (a). From Fig. 2 (b) it can be seen how the correlation coefficient between observed and forecasted water levels goes down as the lead time increased. The 3D hydrodynamic model was built using MIKE 3 FM, a software tool for modelling unsteady three-dimensional flows, which uses a flexible mesh calculation grid taking into account density variations, bathymetry and external forcing, such as meteorology, tidal elevations, currents and other hydrographic conditions ( DHI, 2016 ). Meteorological forcing for the model is obtained from a WRF (Weather Research and Forecasting) limited-area numerical weather prediction model covering Northern Europe with a resolution of 0.1 degrees, which is run by StormGeo in Norway. More information about the calibration, and the model in general, are given in DHI (2011) . A description of the above method for estimating storm surge event time series for climate change analysis is described in Berbel Roman (2014) . Fig. 3 (b) shows the extreme sea water levels which were used as input boundary conditions along the coastline. According to Berbel Roman (2014) , to estimate the expected changes in sea surges in future (up to year 2100), hydrodynamic simulations were carried out, which were driven by the wind and atmospheric pressure results from three Fig. 12. Buildings gradually become visible as the flood water spreads west and northward. The left column shows the plane view of surface topography. The right column shows the corresponding mesh. regional climate models. A general observed extreme event pattern was identified based on the past observed extreme sea-level events. The future extreme water level event time series ( Fig. 3 (b)) were then obtained from the general observed pattern by scaling to a given return period and adding estimates of mean sea level rise and change in storm surge signal. The water level calculated with statistics projections of 100-yr return period considering climate change under a present scenario iden-tified as 'current' and a future scenario identified as '2100 mean' and '2100 upper' depending of the climate factor considered. These extreme sea-level events (current, 2100 mean and 2100 upper) lasted 24 h. As seen in Fig. 3 (b), for a return period of 100 years, the maximum value for the extreme water levels is 1.52 m (current), 2.25 m (2100 mean) and 3.08 m (2100 upper).  Fig. 1 (b)) simulated using MIKE 21 FM and Floodity with a mesh resolution of 20 m, 10 m, and 5 m respectively, based on bathymetric data with/without buildings in a scenario of 100-yr return period rainfall and 2100 upper extreme sea-level event. Extreme precipitation data in this case study is obtained by a frequency analysis, where the rainfall data (over more than 10 years) from 83 stations in Denmark is used. By running a regional statistical extreme model, the intensities of rainfall for the 2-yr, 10-yr and 100-yr return period have been obtained, then multiplied by a factor of 1.1-1.5 due to the impact of climate change in the next 100 years ( Berbel Roman, 2014 ). The designed intensities are shown in Fig. 3 (c). It is assumed that for a given rainfall event, a uniform rainfall is falling over the whole area of the domain (2.2 km × 1.7 km).

Model applications
A series of model simulations using anisotropic-DMO have been carried out to assess the performance of the new flooding model developed here. The extreme seawater level event ( Fig. 3 (b)) which lasts 24 h is used as input boundary condition on the sea boundary ( Fig. 1 (b)), the remaining boundaries are set up as closed (no flow). In these simulations, sea water enters the densely urbanized area from the sea boundary ( Fig. 1 (b)) and the extreme rainfall events ( Fig. 3 (c)) take place simultaneously. Fig. 3 (a) shows the initial water depth within the domain in a scenario of 2100 upper extreme water level. The adapting mesh schemes are listed in Table 1 . The mesh is adapted to ensure an absolute error in the water depth field of 0.01 m and the aspect ratio is 100.  Fig. 4 shows the flood propagation process over the urban area in a scenario of an individual 2100 upper extreme sea-level event. It can be observed that in most of the inundation area, the solutions of water depths obtained from Floodity are in good agreement with those from MIKE 21 FM. The mesh is optimally adapted according to the evolving flow features in time and space, thus providing sufficient mesh resolution where and when it is required (right panel in Fig. 4 ). For example, the fine mesh is located along the flood propagation path while the coarser mesh is used in the areas where inundation has not occurred yet. To further estimate and compare the flood extent, the flood volume during the flood propagation process is calculated ( Fig. 5 ). It is clearly that flood volume obtained from Floodity is a little higher than that from MIKE 21 FM during the whole period, whilst the general trends of both are consistent.

Comparison with MIKE 21 FM results at detector locations
To further evaluate the performance of Floodity using different mesh resolutions, the time series of water depth at detector locations P1, P2, P3 and P4 are plotted in Fig. 6 , in comparison to those from MIKE 21 FM. In Fig. 6 , the blue, green and black lines represent the time series of water depth predicted by the new anisotropic unstructured mesh flooding model Floodity with a minimum mesh size of 20 m, 10 m, and 5 m respectively. It can be seen that a good agreement is achieved between the results from both the fixed and adaptive mesh simulations (except for that with a minimum mesh size of 20 m) at detectors P3 and P4 during the flooding propagation period [6.7, 24] h. However, MIKE 21 FM predicts an earlier flood arrival time at P3 and P4 than that predicted by the Floodity simulations. The results of water depth at detectors P1 and P2 obtained by both the MIKE 21 FM and Floodity simulations are very close to each other when almost the same mesh resolution (10 m) is used over the inundated regions. The detectors P1 and P2 are located within a narrow open channel, where a high resolution mesh (a mesh size smaller than 10 m at least) is required to represent it. We can see that the Floodity simulation with a minimum mesh size of 5 m predicts a deeper water depth at P1 and P2 than that in the simulations with the large mesh size (10 m and 20 m). It proves that more detailed solutions can be obtained in the local areas when using an adaptive mesh instead of a fixed mesh. We also note that using anisotropic-DMO, Floodity is capable of providing reasonable results even with use of a coarser mesh (20 m).
For comparison purposes, MIKE 21 FM results from the model updated from Berbel Roman (2014) are used as a reference solution in this study. The original model domain in Berbel Roman (2014) has a dimension of 2.3 km × 7.5 km, where the northeastern region is selected as our computational domain ( Fig. 3 (a)). Berbel Roman (2014) divided the domain in 9 sub-regions with the element size in range of [10,100] m , where large elements were used to represent the surroundings of harbours, the coastline and the train tracks were represented with smaller elements. Thus, a flexible mesh with the element size in range of [10,100] m is used in the MIKE 21 FM simulations while the adaptive meshes with a minimum mesh size of 20 m, 10 m, and 5 m are used in Floodity. A comparison of water depth results using the adaptive (Floodity) and fixed (MIKE 21 FM) unstructured mesh has been carried out.

Joint flooding events
The newly developed adaptive mesh flooding model is further used for simulating floods under the combined impacts of the events (extreme rainfall and sea level). Fig. 7 shows the results of water depth from MIKE 21 FM and Floodity with a minimum adaptive mesh size of 10 m at time levels = 15 h . It presents the flood propagation process over the urban area in scenarios of an individual 2100 upper extreme sea-level event and a joint event with 100-yr return period rainfall respectively. Due to the effect of rainfall, the joint event has larger flood areas than the individual event, as shown in the areas marked with rectangles in Fig. 7 . Figs. 8 and 9 show the time series of water depth at four detector locations predicted by these models respectively. The total rainfall amount in 24 h is 28.1 mm, 45.7 mm and 82.1 mm within the 2.2 km × 1.7 km study area for the 2-yr, 10-yr, 100-yr return period respectively. Thus, in comparison to the extreme sea level, the rainfalls have a relatively smaller impact on the inundation extent. Again, it can be observed that under the scenarios of joint flooding events, the water depths and velocity obtained from Floodity using the minimum mesh size of 10 m are in good agreement with those from MIKE 21 FM simulations.

Impact of buildings on flooding simulations
The use of anisotropic-DMO in flooding modelling can better capture the evolving flow features and topographic features (buildings, rivers and streets), thus providing improved accurate flooding prediction. To further demonstrate the capability of the flooding model, Floodity has been applied to the joint flooding events with the bathymetric data including buildings, which were represented as impervious obstacles blocking the flow path. Fig. 10 presents the flooding map over the urban area at time level = 15 h in scenarios of an individual and joint flooding event respectively. It can be observed that Floodity results with a mesh resolution 5 m have a larger inundation extent than MIKE 21 FM results and present more details of topographic features, including buildings and channels. The details of roads, buildings and channels can be observed clearly with an increased mesh resolution around them. Fig. 11 provides the details of the areas marked with rectangles in Fig. 10 . It is noted that the information of roads and channels has been lost in the MIKE 21 FM simulations. However, with the use of anisotropic-DMO, Floodity can capture the details of topographic features even with almost the same mesh resolution (10 m) as that used in the MIKE simulation. The bottom in Fig. 11 shows more details of the anisotropic unstructured meshes, where the adapted anisotropic elements are placed under the limitation of the aspect ratio of elements.
In these simulations, the topographical data (digital elevation data) is available with a high resolution of 1.6 m. The availability of highresolution topographical data is important for the accurate numerical simulation of urban food inundation. However, high-resolution topographical data requires a high computational effort, thus, resulting in a computationally demanding flood modelling. Using anisotropic-DMO, the topographical data over the domain is obtained by interpolating the high resolution (1.6 m) data onto the adapted mesh at each time level. Therefore, the high-resolution topographical data is only used in the flooded region while the low-resolution data is used in the rest of the domain, thus reducing the computational cost. In addition, the details of buildings can be represented accurately as the flooding water spreads across the domain (see Figs. 10 and 12 ). Fig. 13 shows the time series of water depth at four detector locations simulated using MIKE 21 FM and Floodity with three different meshes in a scenario of a joint flooding event from extreme sea level and rainfall for a 100-yr return period. It can be observed that due to the impact of buildings, the water depths obtained from Floodity have a slight difference from that of MIKE 21 FM. A deeper water depth at P1 and P2 locations is predicted when using the 5 m resolution adaptive mesh than that using the fixed mesh (MIKE 21 FM), similarly as which is shown in Fig. 6 . This is due to the fact that the detectors P1 and P2 are located at the channel, thus having a larger water depth. Again this proves that the adaptive mesh flooding model can provide relatively accurate predictions. Also note that there is an arrival time lag at detectors if the impact of buildings is considered in the flooding simulations. As a result, the arrival time at P3 and P4 has nearly 140 min' time lag in Floodity with a mesh resolution of 5 m results, in comparison to MIKE 21 FM results.

Sensitivity to the forcing inputs and mesh resolution
Extreme joint flooding is the product of a wide range of interacting processes. Here the uncertainties from the forcing inputs are the extreme sea levels and rainfall. In addition, the mesh resolution is one of critical parameters in flooding modelling. In this section, sensitivity analysis of flood volume over inundated areas to extreme sea levels, rainfall and mesh resolution has been explored and shown in Fig. 3 (a), (b) and (c), respectively.
The impact of the incoming sea levels on flooding results has been investigated and the corresponding results in Fig. 14 (a). There are three scenarios of individual extreme sea-level events shown in Fig. 3 (b), where the maximum value for the extreme incoming water levels is 1.52 m (current), 2.25 m (2100 mean) and 3.08 m (2100 upper) peaking at = 9 . 5 h . In Fig. 14 (a), we can see that the flood volume become large with an increased incoming wave level. The flood volume peaks approximately at = 9 . 5 h when the incoming wave is reaching its extreme. There is a slight time lag in the peak of flood volumes in the scenario of the current extreme sea-level event.
Further investigation of the effect of rainfall on flood volume has been undertaken in the scenario of the joint flood event. A comparison of flood volumes between the scenarios of the individual and joint flood events is provided in Fig. 14 (b). The solid line is the flood volume during the individual extreme sea-level event (2100 upper) while the dashed line represents the joint flood event with the 100-yr return period rainfall event during the rainfall period = 13 − 16 h . The influence of rainfall is reflected by the divergence of the flood volume during rainfall. The largest difference in flood volume is 1.10 × 10 5 m 3 at = 14 . 08 h . Fig. 14 (c) presents the flood volume results from the simulations using different mesh resolutions. It is found that the peak values of flood volume results are very close in all simulations, while the peak time differs greatly when using the mesh resolutions of 20 m (dotted blue line) and 10 m (dotted green line) as well as 5 m (solid black line). The reason for this is that the blocking effect of buildings cannot be represented in flood modelling with use of coarse mesh resolutions (20 m, here) due to the failure of capturing the details of buildings.
In general, the flood volume results are sensitive to both incoming sea levels and rainfall. However, in joint flood events, the effect of rainfall is relatively small in comparison to extreme incoming waves. In flooding modelling the mesh resolution is the key to capture the details of complex topography, for example, buildings and channels. One can see the blocking effect of buildings only when the buildings are captured with high mesh resolutions.

Performance of floodity modelling
As seen in Fig. 15 and Table 2 , an unstructured mesh with 20 m resolution generated by Gmsh ( Geuzaine and Remacle, 2009 ) consists of 13,033 nodes and 25,666 unstructured triangle elements, while 10 m resolution consists of 51,367 nodes and 101,939 elements, and 5 m resolution mesh contains 204,980 nodes and 408,373 elements. These meshes are used during the whole simulation period for fixed unstructured mesh modelling, and used as the initial mesh for the adaptive mesh simulations. After first adapting the mesh, the numbers of nodes and elements used for adaptive mesh of 20 m, 10 m, and 5 m resolution are reduced by 74%, 82% and 88% respectively, then gradually increased during = [1, 600] min and decreased during = [600, 1440] min as the flooding water retreats. Above all, Floodity is less computationally expensive than fixed unstructured mesh modelling. The use of adaptive unstructured meshes improves the computational efficiency. To further reduce the computational cost, various numerical techniques can be adopted in the flood model Floodity, for example, parallel computing using MPI.

Conclusions
Realising the importance of flood coincidence risk assessments, we have further developed the adaptive unstructured mesh flooding model (Floodity, Hu et al., 2018 ) for the joint urban flood events caused by multiple sources (extreme rainfall and sea-level events) and successfully applied to Greve in Denmark. By introducing the anisotropic-DMO technique, the features of flooding flows (local flows around the buildings or the wetting and drying front, for example) are able to be better captured while reducing computational cost without sacrificing accuracy of flooding simulations. With a unique combination of anisotropic-DMO and high-resolution Digital Terrain Model (DTM) data, the complex urban topography can be accurately represented when/where needed by increasing the mesh resolution (around the buildings, for example) dynamically when the flooding water spreads over the urban area. This new Floodity model has been applied to several flooding scenarios that happened in Greve, Denmark, where the flood is induced by different combinations of extreme incoming sea levels and rainfall. A comparison between Floodity and MIKE 21 FM results has been undertaken. It has been found that Floodity is able to provide relatively accurate results while the computational cost is reduced by 20-88% in comparison to fixed mesh models. To assess uncertainties in model predictions, the sensitivity of flood volumes to extreme sea levels and rainfalls has been explored. In joint flood events, we found that the flood volume over the inundated area is more sensitive to sea levels than rainfall. Extreme sealevel events with the higher peak water levels induce higher peak flood volume while the impact of rainfall is relatively small. The sensitivity of flood results to the mesh resolution is also investigated. In flood modelling, the blocking effect of buildings on the peak time of flood volumes can be seen only when using high resolution meshes and Digital Terrain Model data.
Flood modelling is a complex and parametric problem. The input uncertainty is one of the main sources of uncertainty. In this paper, we mainly focused on the simulation of flooding from multiple sources. We have done some basic sensitivity analysis. Given its complexity, in future, we will further carry out uncertainty analysis using advanced numerical techniques, for example, the adjoint sensitivity and uncertainty analysis ( Cacuci et al., 2003 ). In this work, the effect of rainfall is relatively small compared to incoming waves. So infiltration is not taken into account here, namely all amount of rainfall water becomes ponded water on ground surface. We will further introduce infiltration rate in future work. Due to the lack of optimization of codes, the CPU time is not demonstrated here. Instead, we have demonstrated the computational cost is significantly reduced by the decrease of the number of nodes used while the accuracy remains the same or better than that in fixed mesh modelling. In future work we will focus on the optimization of codes (data structures).