Digital Surface Model (DSM) Construction and Flood Hazard Simulation for Development Plans in Naga City, Philippines

Muhammad Zulkarnain Abd Rahman
Department of Remote Sensing
Faculty of Geoinformation Science and Engineering
Universiti Teknologi Malaysia

Dinand Alkema
International Institute for Geoinformation Science and Earth Observation (ITC),
Enschede, The Netherlands

Abstract
Recent development in hydrodynamic modelling requires accurate and detailed terrain model. One of the main problems is frequent changes of land use in major cities, where frequent updating of the digital terrain model (DTM) for flood modelling might be needed. This paper presents an example of assessing the impact of developments on changes in flood behaviours in Naga City, the Philippines. The elevation data is constructed through integrating various elevation data derived from many sources. The development impact assessment is implemented by detailed observation on changes in flood hazard areas.

The geostatistical approach is used to investigate the effect of integrating multi-sources of elevation data by evaluating the nugget values. Four interpolation methods were used, namely Australian National University’s Digital Elevation Model algorithm (ANUDEM), Kriging, Polynomial and Triangulated Irregular Network (TIN). The assessments are based on the Percentile Vertical Accuracy Assessment, error point’s distribution and visual assessment of the DTMs. As a result, the kriging interpolation method has produced the best DTM and it full-filled the requirements for hydrodynamic flood modelling purpose. Two sets of Digital Surface Model (DSM) were constructed to illustrate the situation of Naga City, before and after the developments. The 1D2D SOBEK flood model was used to simulate flood events for 2, 5, 10 and 17.5 years return period flood.

Through this study, it was proved that by simply elevating ground terrain only can solve the flood problem in a particular area. Unfortunately, the flood problem is transferred to another area.

1. INTRODUCTION
Rapid and uncontrolled urbanization in developing countries has become one of the major issues in flood hazard and risk management. This is certainly one of the major environmental problems in the developing world, today and in the years to come. Inevitably, urbanization has great influence on rainfall runoff and flood behaviour. The flow of the floodwater becomes complex due to complicated buildings distribution and structures in an urban area. Excessive water from heavy rainfall easily converted to run-off over paved surface, and due to improper urbanization planning, water will accumulate and increase the potential of flooding. Thus impose great challenge in forecasting and predicting the potential damages. In most cases, improper development in developing countries might ignore above mentioned impact on flood hazard and risk to the surrounding community.

In this case, the 2D flood modeling requires detailed and accurate terrain model, which quality depends on the acquisition techniques and the characteristics of terrain under investigation. The main objective of this research is to generate DTM and DSM of the study area with consideration on current and future developments, followed by simulating the flood events and development impact assessment.

2. STUDY SITE
Naga City is located in Bicol region, at the south-eastern tip of the Philippine island of Luzon. Naga City, located about 377 km to the south of Manila, is well known as a fast-growing area (see figure 2.1). Naga City has the largest population among 35 municipalities in Camarines Sur, which population covers about 8.9 percent of the total population of the province (Naga City Government Philippines Business for Social Progress, 2001). Naga is considered the heart of the Bicol region, consists of 27 Barangays on the land of 7,748 hectares. The main portion of The Naga city is located in low and flat topography that usually inundated by flood when water from the Naga and Bicol River overflow. Thus, substantial discharge and heavy rainfall during monsoon commonly causes severe flood in this city.

Figure 2.1: Naga City

3. DATA COLLECTION RESEARCH METHODOLOGY
In this research the data collection is divided into 4 major groups as follow.

1. Elevation or topographical data
2. Landuse or landcover
3. Recent and future developments
4. Rainfall and floodwater depth during the Super Typhoon Nanmadol.

Additional elevation information were collected to fill gaps in the available data and to update terrain information due to recent and future developments. Other necessary data for instance Landuse or Landcover, Rainfall, flood depth and extent and flood risk perception are needed in flood modeling and development impact assessment. This study is divided into 4 main phases, namely, 1) Data preparation and analysis, 2) DTM and DSM modelling, 3) flood model calibration and modelling and 4) development impact assessment (see figure 3.1). The first phase focuses on elevation data preparation, analysis and integration. The second phase of the research methodology concerns on the construction of DTM and DSM based on different sources of elevation data, which were derived from both primary and secondary sources. The primary data collection aims at filling the gaps in the available data and to update terrain changes in the study area as a result of recent developments. In general, the elevation data is derived in various forms, for instance points, line and polygons, and these data are then aggregated into ground terrain and man-made features. Further aggregation is made to separate the elevation data into two terrain situations; current and future situations. The DTMs are produced using different terrain interpolation methods. The best product is selected and integrated with man –made features to produce DSMs of the study area.

The 1D2D SOBEK flood model is used to simulate 5 recurrence intervals flood events. The flood calibration is made base on flood depth information derived from recent field observations (this data were collected by Saut Sagala and Peters Guarin Graciela) after the flood event caused by the Super typhoon Nanmadol (with an equivalent to 10 years return period flood). The flood calibration is based on two aspects; surface roughness and building structure. The calibrated surface roughness and the suitable building representation will be used for further flood modelling. The surface roughness value of the study is based on landuse or landcover. This updated information is used together with development plans to create recent and future landuse or landcover in Naga City.

The final phase emphasizes on development impact assessments based on detailed investigation of flood characteristics before and after the development. This is supported by additional assessment focuses on changes of flood hazard areas for current and future situation of Naga City. The definition of flood hazard is based on the combination of flood velocity and flood depth.

Figure 3.1: Overall research methodology

3.1 DTM and DSM construction
DTM and DSM of the study area were generated in 4 major steps; 1) elevation data preparation and analysis, 2) elevation data interpolation, 3) accuracy assessment and reporting and finally 4) integrating natural terrain with man made terrain to produce DSM.

3.1.1 Elevation data preparation and analysis
The elevation dataset for DTM generation is derived through the integration of various elevation data sources and these data vary in both horizontal and vertical accuracies (see appendix 1). The arrangement of elevation data will influence the shape of the variogram model. In this case, the clustered elevation data, which were found around dunes, road embankments and other local abrupt changes in topography, were removed from the dataset (Blomgren, 1999). Wilson and Atkinson (2003) in their research, “Prediction the uncertainty of DEM on flood inundation modelling”, used the ordinary kriging to interpolate the elevation data in the floodplain area. The elevation data was the combination between the contour lines and the elevation points that were derived from GPS measurement. The original experimental semi-variogram of the contour lines had quite general shape. However this general shape or trend was reduced (increased variance at shorter lags than globally) when the elevation points derived from GPS measurements were added to the dataset.

Ten set elevation data were used and integrated to produce DTM. These datasets vary in scale and contour interval which remarks difference in horizontal (planimetric) and vertical accuracies. The problem of integrating elevation data from different sources with different scales and accuracies lies on the fact that the elevation values in the combined dataset may lie close to each other. The challenge is to identify an appropriate approach to prioritize the datasets, to identify which of those datasets represent the true terrain elevation and to combine the entire datasets. Thus, the datasets are prioritized with 2 steps; 1) Prioritization based on Nominal horizontal and vertical accuracies (based on the National Standard Data Accuracy (NSSDA)), 2) Prioritization based on data forms (spot heights and contour lines) and production date (see table 3.1).

Table 3.1: Available elevation data sources for DTM generation

S.No. Source of data Priority level based on the nominal accuracy Priority level based on the 2nd selection step Final priority level

The contour lines are converted to points and then combined with other point form data (spot heights). Data with higher priority score would replace the lower priority score data. The replacement is done when 2 or more elevation points fall within 3 meters radius. The effect of the integration of the multi-sources elevation data is assessed by means of semi-variogram analysis. The assumption is points that are close together should have less difference or high autocorrelation. Thus, high nugget value would reveal strong effect of disagreement between the elevation datasets. Certainly, the nugget effect could also attribute to the complexity of terrain features. However it was found that, the effect of data disagreement still appear when datasets with complex geomorphological features were removed. Several attempts with different integration method were used (see table 3.2) to decrease the nugget value. However, for the sake of the terrain complexity information, the nugget value was reduced from 2.9 to 2.2.

Table 3.2: The value of Semi-variogram model parameters for each dataset; the 2nd and the 3rd datasets will be used for the DTM interpolation

DatasetNuggetSillRangeLagModelNumber of points

At this stage, the integrated elevation datasets with low nugget value are assumed to have less disagreement between elevation dataset, less overlapping dataset, good elevation data in representing the real terrain and inevitably contain degraded complex terrain features.

3.1.2 Elevation data interpolation
The DTM interpolation method should in general preserve the detailed terrain information while reducing the effect datasets disagreement. The interpolation of the ground terrain is done with 4 interpolation techniques; 1) Kriging, 2) TIN,3) Polynomial trend surface and 4) ANUDEM (see figure 3.2).

A Geostatistical interpolation or kriging interpolation method is similar to a probabilistic interpolation technique, in which the weights are derived from the surrounding sample points. However, the weights are not only based on the distance, but also on the strength of the overall correlation among the measured points (Maune, 2001). The basic interpolation assumption is that, values at a short distance are more likely to be similar than at a larger distance. Demirhan et al., (2003) in his study had focused on the performance of several interpolation methods with presence of noise and sampling pattern. He had pointed out that the “Ordinary Kriging is the most robust interpolations method against noise”. On the other hand, the deterministic DTM interpolation approach tries to fit a mathematical function to a set of elevation samples of known coordinate (x and y). This interpolation technique could be done either through an exact interpolation or smooth interpolation (Meijerink et al., 1994).

Figure 3.2: Natural terrain derived from 6th Polynomial degree (a), TIN-based terrain modelling (b), ANUDEM (c), Ordinary Kriging (100 m conversion interval) (d), and Ordinary Kriging (5 m average block) (e), visualized in 3D

The smooth interpolation method is suitable with the assumption that point measurement data are regarded as true value or with less error. On the other hand, with coarse data accuracy, the smooth interpolation scheme might be the best way to level out the error to some degree. The fitting process uses various types of mathematical functions usually known as polynomial functions at a certain degree of complexity to fit the surface through the sample points.

The nature of TIN modelling is splitting the surface into triangular element planes (Meijerink et al., 1994). More detailed definition, “TIN is a digital terrain model based on irregular array of points which forms a sheet of non-overlapping contiguous triangle facets” (Maune, 2001). The next interpolation technique is ANUDEM. ANUDEM is a software package known as the Australian National University Digital Elevation Model developed by Hutchinson (Geodata and Geoscience Australia, 2002). It was designed and optimised to create a hydrologically correct terrain models. “It’s unique in both input and output for building a good terrain model”(Maune, 2001). The input data is not only confined to point data, but also lines which represent streams and ridges for drainage, and polygons as a lake boundary to produce a DEM that is virtually free of spurious sinks and pits

3.1.3 DTM quality assessment and reporting
The DTM quality assessment encompasses 3 main parts; 1) DTM accuracy assessment 2) DTM errors distribution and 3) General gemorphological of the study area. The DTM accuracy assessment uses Percentile Accuracy Assessment method introduced by Maune (2001) which is suitable for non-normally distributed residual data. According to table 3.3, ANUDEM interpolation method has produced the most accurate DTM with the maximum error 0.98 m at 80 percent of the total residual. The DTM produced by the Ordinary Kriging with 5 m block average has the second lowest error with maximum error 1.00 m. followed by DTM produced by Ordinary krigging on 100 m contour lines to point conversion and finally the polynomial interpolation methods.

The flood modelling gives more focus in the low-lying areas, where most of the commercial areas are located. Thus, the next assessment focuses on the number of point with errors more than 1.0 within area of interest (see appendix 2). It was found that, TIN-based surface modelling has the least number of points with error more than 1.0 m and followed by the Ordinary Kriging with 5 m block average, the Ordinary Kriging with 100 m contour lines to points’ conversion interval, ANUDEM and Polynomial. The final stage in DTM quality assessment focuses on the capability of DTMs in representing the real landscape of the study area. According to figure 3.2, the DTM produced by the Ordinary Kriging with 5 meters elevation block average is able to represent the real landscape of the study area. The DTM has shown in general, back-swamp areas, natural levee and elevated areas for bridge construction. However, another DTM produced for instance by ANUDEM and TIN contain large number of artificial pits. According to Raaflaub and Collins (2005) artificial pits or sink features in DTM are hydrologically serious problem. On the other hand, DTMs produced by the Polynomial interpolation method shows smooth surface and failed to represent the general landscape of the study area. Therefore based on above judgements, the DTM produced by the Ordinary Kriging with 5 meter elevation block average is used for DSM construction.

Table 3.3: Summary on DTMs quality assessment

Interpolation methodAccuracy assessmentNumber of errors with more than 1.0 m in the area of interestVisual assessment

3.1.4 Integrating natural terrain and man-made terrain
The DSM of the study area is constructed through combining the simulated man-made terrain and DTM. The DSMs are used to represent current and future situations of Naga City (see figure 3.3). The man-made features are simulated by integrating spatial data (e.g. roads, drainages, landuse or landcover), building footprints and detailed development plans.

Figure 3.3: DSM represents current situation (a) and future situation of Naga City (b)

3.2 Flood modeling and flood model calibration
The 1D2D SOBEK flood model is used to simulate flood behavior which dominantly occurred due to substantial discharge of the Naga and thr Bicol River. This model requires incorporation of detailed floodplain terrain model and river cross sections, which both parts are represented by 2D and 1D model. The flood modeling covers 4 flood recurrence intervals namely 2, 5, 10 and 17.5 years return period flood. Flood depth and extent information of flood event triggered by the Supertyphoon Nanmadol is used in flood model calibration. Previous study had confirmed that the magnitude of this flood event is comparable to 10 years return period flood. The flood model calibration focuses on calibrating 2 parameters, 2 sets of building structures and 3 sets of surface roughness. In the flood model calibration process, it was found that, water velocity of rough surface building structure is higher compared to the solid block. On the other hand, the rough surface building structures allows water to flow much farther and increase the flood extent. The obstruction has increased the flood depth and subsequently increase flood velocity at the edge of the buildings. Besides, there is no significant difference in flood depth for flood models with different set of manning coefficients. Buildings in Naga City are quite dense for both residential and commercial areas, thus solid-block building structures is selected for further flood modelling.

3.3 Development impact assessment
Development impact assessment makes use of changes in flood hazard areas in delivering messages on the changes of flood impact on current and future situations of Naga City.

3.3.1 Changes in flood hazard areas
Flood hazard can be expressed by various combination of flood characteristics, for instance, flood depth and flood velocity based on research by Ramsbottom et al. (2003) (see figure 3.4). Hazard here specifically refers to the wading hazard of adults, children and also selected vehicles during flood. In this assessment, flood hazard maps for each simulation hour are created and integrated based on the worst case of hazard level at each pixel (see figure 3.5).

Figure 3.4: Graph that combines flood velocity and flood depth to define level of flood hazard (Ramsbottom et al., 2003) (a) and Major developments in Naga City (b)

Figure 3.5: Flood hazard distribution for current (left) and future (right) situations for 17.5 years, 10 years and 5 years return period flood in Naga City

5.0 CONCLUSION
Step by step multi-sources elevation data integration is one of the cheapest approaches in constructing terrain model without precise and expensive elevation data acquisition techniques for instance, LiDAR, Radar Interferometry and Aerial Photos. However, the integration of all the available elevation datasets should be made very carefully. In this case, the terrain updating process can be made through the conventional levelling and the compilation of recent and future development plans to be added onto the natural ground terrain model.

Figure 3.6 shows detailed changes in flood hazard areas in Naga City for 17.5 and 10 years return period flood. According to figure 3.6 (b), for 10 years return period flood, it is about 92 percent of the areas with the positive impact are located in commercial area. Meanwhile, areas with negative impact are mostly in the residential and other commercial areas (figure 3.6 (a)). On the other hand, for 17.5 years return period flood (figure 3.6 (d)), commercial areas account about 96 percent of the total areas with positive impact. Besides, it is about 34 percent of the agriculture areas, 21 percent of the residential areas and 13.88 percent of the commercial areas are classified in negative impact areas.

It must be realised that the development impact assessment doesn’t intend to oppose any developments in Naga City. This simulation results would give better ideas and understanding on flood behaviour which subsequently aids flood management and mitigation processes. Development is necessary to provide current necessities and needs of the future generations. However, sustainable development may not be achieved without any effort in reducing the intensity of hazard as a result of developments.

Figure 3.6: Areas classified into negative and positive impact (based on the hazard map) for 10 years return period flood (a) and (b) and 17.5 years return period flood (c) and (d)

7.0 Appendices
Appendix 1: Source of elevation data for DTM generation

Appendix 2: Errors distribution for TIN (a), Ordinary Kriging with 5 m block average (b), Ordinary Kriging (100 m line to point conversion) (c) and ANUDEM (d); TIN has the least error within the study area, followed by block ordinary kriging with 5 m, Ordinary Krigging for 50 m line to point conversion interval and ANUDEM

References

• Blomgren, S. 1999. A digital elevation model for estimating flooding scenarios at the Falsterbo Peninsula. Environmental Modelling & Software 14:579-587.
• Demirhan, M., A. Ozpinar, and L. Ozdamar. 2003. Performance evaluation of spatial interpolation methods in presence of noise. International Journal of Remote Sensing 24:1237-1258.
• Geo-computation (ed.) 2003. 7th International Conference on GeoComputation, University of Southampton, United Kingdom. 8 – 10 September 2003. GeoComputation.
• Geodata and Geoscience Australia. 2002. Data Use Guide; A digital elevation model of Australia with a grid spacing of nine seconds in Longitude and Latitute [Online] https://chapters.marssociety.org/canada/
  expedition-mars.org/ExpeditionTwo/files/maps/DEM_userguide.pdf (verified 22/11/2005).
• Maune, D.F. 2001. Digital elevation model technologies and applications : the DEM users manual The American Society for Photogrammetry and Remote Sensing (ASPRS), Bethesda.
• Meijerink, A.M.J., J.A.M. de Brouwer, C.M. Mannaerts, and C.R. Valenzuela. 1994. Introduction to the use of geographic information systems for practical hydrology : IHP-IV M 2.3 ITC, Enschede.
• Naga City Government Philippines Business for Social Progress. 2001. Naga City Disaster Mitigation Plan. Working paper. Asian Disaster Preparedness Center, Bangkok.
• Raaflaub, L.D., and M.J. Collins. 2005. The effect of error in gridded digital elevation models on the estimation of topographic parameters. Environmental Modelling & Software.
• Ramsbottom, D., P. Floyd, and E. Penning-Rowsell. 2003. Defra / Environment Agency, Flood and Coastal Defence R&D Programme. Defra – Flood Management Division, London.
• Smith, K. 2001. Environmental hazards : assessing risk and reducing disaster. Third edition ed. Routledge, London.