Influence of elevation modelling on hydrodynamic simulations of a tidally-dominated estuary

Journal of Hydrology 497 (2013) 152–164

Contents lists available at SciVerse ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Influence of elevation modelling on hydrodynamic simulations of a tidally-dominated estuary Ana Paula Falcão a,b,⇑, Andrea Mazzolari a,c, Alexandre B. Gonçalves a,b, Maria Amélia V.C. Araújo a,c, António Trigo-Teixeira a,c a b c

Instituto Superior Técnico, Technical University of Lisbon, Av. Rovisco Pais, Lisboa, Portugal ICIST, Instituto de Engenharia de Estruturas, Território e Construção, Lisboa, Portugal CEHIDRO, Centro de Estudos de Hidrossistemas, Lisboa, Portugal

a r t i c l e

i n f o

Article history: Received 18 August 2012 Received in revised form 6 May 2013 Accepted 26 May 2013 Available online 4 June 2013 This manuscript was handled by Andras Bardossy, Editor-in-Chief, with the assistance of Uwe Haberlandt, Associate Editor Keywords: Spatial interpolation Digital elevation model resolution Hydrodynamic simulation in estuary Satellite positioning High resolution geoid undulation model Wet and dry

s u m m a r y Hydrodynamic simulation of estuaries requires a single digital elevation model (DEM) resulting from merging of both topographic and bathymetric data. These two datasets are usually produced using different technologies, co-ordinate systems and datums. Intertidal data in particular are often lacking due to the difficulty of data acquisition using conventional survey techniques. This paper presents a fast, accurate and low-cost methodology to fill this gap and highlights the effect of the digital elevation model characteristics, such as the interpolation method and spatial resolution, on modelled water levels and flooded areas. The Lima river estuary, located in North-western Portugal, is used as a case study. Validation tests for commonly available spatial interpolators showed ordinary kriging to be the most adequate interpolator. Digital elevation models with regular grids of 5 m and 50 m resolution were used, together with the original (not interpolated) elevation dataset, as input to a finite element hydrodynamic model for astronomic tide simulation. Results indicate that the larger differences between using different elevation models occur at low tide during spring tide, marginally impacting the flood modelling. The effect of a vertical offset of the chart datum with respect to a part of the digital elevation model was finally investigated, showing a limited influence in the determination of the water levels. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction Floods in estuarine areas occur when the normal tidal level is exceeded. These are natural phenomena which are mainly linked to the surge of sea level caused by meteorological factors, such as low atmospheric pressure, high freshwater flow rates or the unusual entrance of water in the estuary due to high run-off events (McInnes et al., 2009). The study of flood events, including the estimation of the highest water level and the flooded area, has to rely on a model that has been previously validated for ordinary flow conditions: in the case of an estuary domain, these are represented by the astronomical tide forcing and the river flow. Simulations made by hydrodynamic models have limitations in considering physical and hydrodynamic conditions, such as the maximum number of elements used to describe the channel geometry and roughness (Merwade et al., 2008a), and these became more evident in estuarine areas, where a series of additional circumstances influencing the flooding process have to be considered ⇑ Corresponding author. Address: Department of Civil Engineering, Architecture and Georesources, Instituto Superior Técnico, Av. Rovisco Pais, 1049-001 Lisboa, Portugal. Tel.: +351 218418328. E-mail address: [email protected] (A.P. Falcão). 0022-1694/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jhydrol.2013.05.045

(shallow waters and intertidal areas, distinct flow directions and complexity of the drainage system). 2D/3D hydrodynamic models require an integrated continuous surface combining elevation from the river bed (bathymetry) with the topography from adjacent areas covering the floodplain. This may eventually be complemented by the inclusion of more detailed data representing features that influence the spread of flooding, such as dykes or roads (Poulter and Halpin, 2008). Results of the hydrodynamic modelling processes are known to be dependent on the properties of the elevation model (Ali et al., 2009; Weaver and Slinn, 2010; Cea and French, 2012). The channel bathymetry, namely the horizontal resolution, is known to affect the results of models, such as the rate, extent and timing in inundation mapping (Hardy et al., 1999; Omer et al., 2003; Horritt et al., 2006; Buttner, 2007; Raber et al., 2007; Poulter and Halpin, 2008; Merwade et al., 2008b). Despite these known effects, changes of scale, which result in distinct horizontal resolutions for the input elevation data, may be necessary to obtain files that do not overload the hydrodynamic model. To create the digital elevation model (DEM) for the simulation stages, data from several sources, formats and acquisition and processing techniques, have to be combined: discrete bathymetry

A.P. Falcão et al. / Journal of Hydrology 497 (2013) 152–164

points, surveyed cross-sections, satellite and aerial elevation models and/or imagery, DEMs obtained from paper or digitalized topographic maps (contours and height spots), topographic features such as thalwegs and artificial constructions such as dykes, piers or culverts. These spatial data may be acquired through a range of techniques, such as aerophotogrammetry (Yamano, 2007; DiGruttolo and Mohamed, 2011), LiDAR (Coveney and Stewart Fotheringham, 2011; Pe’eri and Long, 2011) or videogrammetry (Long, 2005), but the application of those techniques, which provide high-resolution data, is only justified in large areas. All these data have to be integrated and many aspects, related with scale, spatial reference, or the choice of interpolators, have to be addressed. Several studies have identified the choice of interpolation method and spatial resolution as decisive factors on the accuracy of the results (Aguilar et al., 2006; Chaplot et al., 2006; Legleiter and Kyriakidis, 2008; Guo et al., 2010) but rarely applied to the hydrodynamic modelling with data acquired from very different sources, such as positioning via Global Navigation Satellite Systems (GNSS), echo-sounding and topographic techniques. Merwade et al., 2008a present a varied set of GIS techniques for merging the datasets of bathymetry and topography from surrounding areas, and acquired by different techniques, in order to create a continuous river and floodplain elevation model. For this, a series of procedures for interpolating and merging those varied datasets is described, but only focused on DEM production. In one of the few modelling applications to real domains, Cook and Merwade (2009) combine datasets obtained from different river topography sources and 1D/2D model approaches, in order to quantify the differences in the inundation maps. Cea and French (2012) present another application of a hydrodynamic model to an estuary using composite data, although the emphasis is the effect of realistic uncertainty ranges in the various topography/bathymetry zones on the modelled elevations and velocity field.

2. Objectives and methodology The overarching goal of this paper is to evaluate and discuss the influence of elevation modelling on hydrodynamic simulations, focusing on a tidal estuary. To address this, the following specific objectives are considered: firstly, a proposal of an efficient, lowcost and operational method that enables the acquisition of data in the intertidal area is made, in order to build a single elevation model, that integrates the bathymetry and topography data, as required in hydrodynamic simulations; a second goal is to compare the performance of commonly available interpolators to produce such an elevation model; a third purpose is to compare and discuss the results (flooded area and water level) of the ADCIRC hydrodynamic model (Luettich et al., 1992) in an estuarine area when the horizontal resolution of the DEM input to the finite element model changes; finally, it is proposed to analyse the relation between the possible errors in the vertical referencing of the bathymetry and their effect on the computed water levels and inundated areas on the estuary. For the first purpose we follow a methodology that combines bathymetric data with topographic data and GNSS-based points acquired for the intertidal area and for control points near the estuary. Considering the nature of the data collected by the GNSS technique, it is necessary to create a local geoid undulation model. This requires a choice of the spatial interpolation method and a thorough analysis of the results, as the interpolated surfaces can be very diverse. It is also necessary to evaluate the influence of the horizontal resolution in the results of the hydrodynamic modelling, by comparing the effects of distinct DEM regular grids. The same

153

applies for the analysis of the influence of errors in the vertical referencing of the bathymetry.

3. Application 3.1. Description of the area and data The area studied in this paper is the Lima River estuary, in the region of Viana do Castelo, North-western Portugal (Fig. 1). The Lima river estuary has already been subject of different numerical applications, describing tide propagation (Rebordão and Trigo-Teixeira, 2009), salt water intrusion (Pinho and Vieira, 2007) and the effect of wave set-up during a storm surge (Araújo et al., 2011). None of these studies, however, considers the sensitivity of the modelled hydrodynamics to DEM quality and resolution. The studied area extends for approximately 10 km from the Lima river mouth. The river flows about 65 km in Portugal and the city of Viana do Castelo lies in its estuary. The north and south banks are connected by the Eiffel and the highway bridges. In this region, agriculture plays an important role in economy, so the management of the low-lying lands is of utmost importance. The estuary is exposed to semidiurnal tides, which represent the main hydrodynamic forcing action, with a range varying from 1.1 m on neap tides to 3.7 m on spring tides. The river regime is controlled by the presence of two hydroelectric power plants upstream, which maintain an environmental flow throughout the year. The domain can be divided into three main parts: In the first part, closer to the river mouth, the average width of the main river channel is about 300 m. The north bank is fully consolidated by ripraps and revetments, with the presence of harbour docks, a shipyards and a marina. The south bank has some mooring platforms, docks, mixed with farming (dimensions about 300 m); in the middle part, between the two bridges and whose length is about 2500 m, the average width of the channel is about 1000 m and intertidal area is about 500 m in north bank and 300 m in south bank; in the last part, the channel dimension decreases to 300 m and intertidal area has an average width of 100 m on the north bank and on the south bank has an average of 400 m. The intertidal area has a very gentle slope of less than 2%. The implementation of the hydrodynamic model considers two spatial scales, one oceanic and one estuarine (Fig. 5). At the oceanic scale, the bathymetry was retrieved from the Global seafloor database from the Institute of Geography and Planetary Physics (Smith and Sandwell, 1997), integrated with data from the 1:150,000 and 1:1,000,000 Portuguese nautical charts for depths between 4 m and 4000 m. At the estuarine scale, existing datasets (topography and bathymetry) were used to provide a single digital elevation model, but at this scale it was still necessary to acquire additional data. There is a variety of positioning techniques available for data acquisition, from a cost-benefit perspective the choice should take into account factors such as the quality, the size of area and the time spent on their acquisition. For this reason, and by the positive balance between accuracy and cost, it was decided to collect the missing data by GNSS surveying. The topographic dataset comprises contour data and height spots, provided by the Portuguese Army Geographical Institute, from 1:25,000 topographic maps (10 m interval). The bathymetric data were acquired during a surveying campaign with different equipment and procedures: a bathymetry dataset acquired with multibeam sonar (99,399 points), in regular grid format, with 5 m spacing (named Estuary in the paper) and 65,535 point set,

154

A.P. Falcão et al. / Journal of Hydrology 497 (2013) 152–164

Fig. 1. Lima river estuarine area.

Fig. 2. Topographic (contours and height spots) and bathymetric (Estuary and Lima) datasets and acquired GNSS points.

surveyed upstream the Eiffel bridge, with 25 m  10 m spacing, named Lima (see Fig. 2). The height reference for all topographic datasets in Portugal mainland is the tide gauge mean sea level at Cascais, defined by an average of the 1882–1938 tide gauge records. Bathymetric data are relative to the hydrographic zero level (chart datum), defined as 2.0 m below the height reference, as presented in the Portuguese Hydrographic Institute tide tables (Instituto Hidrográfico, n/d).

3.2. Data collection 3.2.1. GNSS data A Leica GS-20 GNSS receiver and a Leica AT501 mono-frequency antenna were used to record 360 positions per point with an interval of 5 s, using the Global Positioning System. A total of 182 points in the intertidal area (Fig. 2) where data were missing were collected in September 8th–10th, 2010. Considering the objectives

155

A.P. Falcão et al. / Journal of Hydrology 497 (2013) 152–164 Table 1 Statistics for positional uncertainty of GNSS surveyed points. Statistic

Planimetry (cm)

Altimetry (cm)

Range Average Standard deviation

14.5 2.5 1.9

20.5 3.3 2.8

of the work and the unavailability of real-time positioning equipment, the final coordinates were obtained with a post-processing procedure, using a local geodetic mark as reference (a short base, about 2 km). Table 1 lists the positional quality values obtained after the post-processing procedure. As the collected coordinates are relative to the WGS84 ellipsoid, it was necessary to convert those to orthometric heights using a local model of geoid undulation, according to:

hðPÞ ffi HðPÞ þ NðPÞ

ð1Þ

where H(P) denotes the orthometric height, h(P) the ellipsoidal height and N(P) the geoid undulation relative to a generic point P. Hence, as N(P) is unknown, there was the need to set up a local geoid undulation model. 3.2.2. High resolution geoid undulation model construction To support the construction of the local geoid undulation model, a sample of three trigonometric marks and seven levelling marks were considered. The local geoid undulation model was obtained by spatial prediction (Fotopoulos, 2005; Nahavandchi and Soltanpour, 2006; Zaletnyik et al., 2007; Falcão et al., 2008). The orthometric heights of this dataset were provided by the National Institute of Cartography (Instituto Geográfico Português – IGP), and the ellipsoidal heights were collected by GNSS equipment – a double-frequency receiver (5 s interval, 720 positions per point, reading GPS and GLONASS frequencies) and a choke ring antenna. The ellipsoidal coordinates were obtained with a postprocessing procedure. The final geoid undulation values, for this dataset, were obtained using expression (1). A polynomial data trend was identified using overfitting criteria (Fotopoulos, 2005) and subtracted to the original values. Residual values were used to predict geoid undulation. According to the characteristics of the residual values (intrinsic stationarity, constant mean and isotropic behaviour), the ordinary kriging technique was used. The exponential covariance model was selected as a theoretical model, given as (Wackernagel, 2003):

CðqÞ ¼ r2 expð3q=aÞ

ð2Þ 2

where q represent the distance between each pair of points, r represents the variance and a represents the effective range. Both these parameters were estimated from empirical variograms. Predicted values (z⁄) were obtained by a weighted average of the observations (z) (Wackernagel, 2003; Reguzzoni et al., 2005; Qiong and Guodong, 2011), according to:

z ¼

n X

ki zi

weighted with elliptical radius (IDWe), thin plate spline with tension (TPS), completely regularized spline (RS) and ordinary kriging (OK) – are a subset of the group also analysed by Merwade et al. (2006) for cross-sections bathymetry. All are implemented in the ArcGIS Desktop 10 software with the Spatial Analyst and the Geostatistical Analyst extensions. On this collection only the IDWe and OK interpolators consider the anisotropic nature of the data. Details and discussions of the various interpolators and of its applications in the construction of bathymetry models can be found in Chaplot et al. (2006), Merwade et al. (2006), Aguilar et al. (2006), Merwade et al., 2008a and Merwade (2009). IDWc and IDWe were parameterised with the most common power value (2) and the 12 nearest points for a maximum radius of 800 m. In the cases of the bathymetry data this value far exceeds the maximum distance between ‘‘neighbour’’ points, but for the topography this valued seemed adequate to manage the larger distances observable in the dataset. The larger semi-axis of the ellipses used in IDWe was aligned with the flow direction, with an anisotropic ratio of 26.6, to consider the greater variability in the transverse direction. TPS and RS were used with the parameters that matched those of the run for IDWc. For OK a first order polynomial detrending was applied to the original dataset after the overfitting criteria analysis, and the residual values were used in the prediction. The exponential model was chosen as the theoretical model and parameters were extracted from the empirical variogram, similarly to the local geoid undulation model estimation procedure. In order to proceed with the validation, the original dataset with 115,597 points combining bathymetry and topography was subdivided into two datasets: a training set containing 80% of all points and a test set with the remaining 20%. As the datasets were large and well distributed spatially, a single random subdivision was made. Two validations were made for each of the interpolators: cross validation and independent validation. Statistical results are presented in Table 2. From these results it is noticeable that both validation procedures present similar statistical results, which was expected due to the large number of input points and the circumstance of being a well distributed dataset. While the mean value of the error for all the five interpolators is small, the range of variance values suggests that OK had the best performance. Hence, it was the selected operator to set up the final combined elevation model. 3.2.4. Final digital elevation model The final digital elevation model was obtained by combining all datasets (contour data and height spots, Estuary, Lima and the GNSS points in the intertidal areas) in three different situations: with linear interpolation of the raw scatter sets, and with spatial prediction, both converted to a regular 5 m spaced grid and to a regular 50 m spaced grid. The data trend was previously modelled by a bilinear function and its values were subtracted from the original values. Ordinary kriging was then applied to this difference map and the final digital elevation model was obtained (Fig. 3).

ð3Þ

i¼1

Values for the weights ki were chosen in order to minimise the prediction error variance. 3.2.3. Comparison of interpolation methods A set of commonly used interpolators was selected to analyse their performance in representing the bathymetry and topography data into a combined dataset. The used interpolators – inverse distance weighted with circular radius (IDWc), inverse distance

Table 2 Statistical results of cross and independent validations for each interpolator. Interpolator

IDWc IDWe TPS RS OK

Cross validation

Independent validation

Mean (m)

St. dev. (m)

Mean (m)

St. dev. (m)

0.06 0.01 0.01 0.03 0.01

1.85 4.27 1.67 1.59 1.32

0.06 0.02 0.02 0.04 0.01

1.99 4.35 1.83 2.02 1.50

156

A.P. Falcão et al. / Journal of Hydrology 497 (2013) 152–164

To test this last hypothesis, and to understand how possible errors in the vertical referencing of the bathymetry have an impact on the simulation results, a practical test was made by introducing an arbitrary vertical difference of +0.25 m to the original elevations of the Lima dataset (see Section 3.3.4). 3.3. Hydrodynamic modelling

Fig. 3. Lima estuary. Final combined digital elevation model.

3.2.5. Validation of elevation The described GNSS-based operational methodology may also be applied to validate the bathymetry datasets. A comparison of the acquired echo-sounding bathymetry and the GPS heights converted onto a common vertical reference was made. A sample of randomly selected 19 points in the north bank (Fig. 4), for which ellipsoidal heights were collected, was used. Results of the comparison with the surrounding bathymetric dataset vary from a minimum of 0.46 m to a maximum of 0.59 m, with a mean value of 0.20 m and standard deviation of 0.25 m. The main explanation for this difference is related with the distinct dates of acquisition, as bathymetry collection and field work datasets have an interval of four years (2006–2010). A minor effect that might contribute for these differences is related with the swampy terrain, on which the operator with the antenna pole stood during the minutes of data collection. This caused the antenna pole to sink up to 5 cm until harder soil is found underneath, as no broad-based survey pole was available. Another explanation may be related with the existence of a systematic error in the bathymetry dataset.

3.3.1. Model description ADCIRC (Luettich et al., 1992; Luettich and Westerink, 2004) is a two dimensional finite element model widely applied to ocean and coastal domains. ADCIRC solves the depth integrated barotropic shallow water equations under the hypothesis of incompressibility and hydrostatic pressure. The model has been shown to provide accurate solutions for domains encompassing different spatial and hydrodynamic scales (Blain et al., 1994; Fortunato et al., 1999; Dietrich et al., 2011). ADCIRC includes a wet and dry algorithm (Luettich and Westerink, 1999; Blain et al., 2010) for reproducing the covering and uncovering of inundated areas, based on an element-by-element evaluation of the water elevation and velocity. In the present modelling study, the computational domain is divided into two spatial scales, one oceanic and one local. The model is first run for the oceanic domain extending for around 615 km westward of the Iberian Peninsula (Fig. 5a). The water elevations calculated in the nearshore in front of the Lima river mouth are imposed as ocean boundary condition for the estuary domain (Fig. 3b). The estuary boundary is approximately located along the +4 m isoline above mean sea level (MSL), in order to cover the maximum flooded area at high tide, being the maximum high water spring +1.95 m above MSL. The same figure shows the deployment sites of 3 Nortek 2 MHz AquaDopp Current Profilers

Fig. 4. Validation point set (in magenta) and Lima dataset (in blue). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

A.P. Falcão et al. / Journal of Hydrology 497 (2013) 152–164

157

Fig. 5. (a) Ocean domain and (b) Lima estuary domain, North of Portugal (depths are referred to the mean sea level).

(ADCPs), which registered the water elevation and the depth average current time series used in the calibration of the estuary model. The water elevation time series was corrected for variations in atmospheric pressure. The domain extends upstream until the location of ADCP 3. As a first approximation, no river discharges are introduced in the simulations. For a good model performance, two crucial requirements have to be considered: the definition of a sufficiently accurate DEM (Cea and French, 2012), and the design of a finite element mesh which is able to describe the complex bottom morphology: a proper mesh resolution should be guaranteed in the near-shore region (Jones and Davies, 2007) and along hydraulic and vegetation features likely to influence the flood propagation path (Cobby et al., 2003; Feyen et al., 2006; Westerink et al., 2008). Only when these requirements are met, the model can be validated for the reproduction of extreme events. Two different finite element meshes discretizing the domain are tested: M5, using the 5 m resolution DEM grid, and M50, using the 50 m resolution DEM grid. In addition, a third test is considered, where the raw dataset of the estuary bathymetry, obtained by merging the Estuary, Lima and contour data scattersets, is linearly interpolated, without kriging treatment, to the finite element mesh M5. For this last test, the GNSS survey points have not been introduced in the resulting DEM. Some details of the meshing process and of the model validation tests are reported after.

3.3.2. Mesh generation The multi-criteria meshing strategy (Mazzolari et al., 2012) is adopted here: a series of a priori criteria, which may come from physical, hydrodynamic, geometrical or other input variables, is defined over the domain. The target node distribution of the prospective mesh is determined by the spatial overlay of all the chosen criteria in a piecewise linear Node Spacing Function (NSF). An advancing front algorithm eventually performs the domain discretization and creates the triangular finite element mesh to be used inside ADCIRC. For the ocean domain, a function which keeps constant the ratio between tide wavelength and mesh size (Le Provost and Vincent, 1986) is used. For the estuary domain the inverse of the bathymetry gradient and the horizontal length scale are considered simultaneously. The first one is expressed as:

Dx1 ðx; yÞ ¼

a ~hðx; yÞj jr

ð4Þ

~hðx; yÞ is the bathymetry gradient and a a coefficient of where r proportionality. Eq. (4) expresses the increase in the mesh resolution in correspondence of rapidly varying bottom slopes. The horizontal length scale criterion, referred as Dx2, sets the minimum element size imposed along the cross section of the principal morphological elements of the estuary, like the main and secondary river channels, the harbor navigation channel and the visible emerged bed forms. These values are mainly empirical, derived from the

158

A.P. Falcão et al. / Journal of Hydrology 497 (2013) 152–164

Table 3 Element size values of Dx2 adopted for M5 and M50. Estuary element

Dx2 (m) applied to M5

Dx2 (m) applied to M50

Harbor navigation channel Main river channel Secondary river channel River sand and gravel bars Floodplains

20

40

20–25 30–35 30–60 >30

40 50–55 50–65 >45

available bathymetry resolution and from the modeller’s experience of similar meshing studies. The Dx2 values used for M5 and for M50 are reported in Table 3. The final NSF is defined as the minimum value between Dx1 and Dx2. M5, obtained by applying the mentioned criteria and interpolating from the 5 m DEM grid, has 17916 nodes and 34762 elements, with an enhanced resolution of the main river and navigation channel (Fig. 6). M50, which applies the same criteria to the 50 m DEM grid (Fig. 7), has 8229 nodes and 15667 elements. 3.3.3. Model calibration The model calibration is composed by an ocean simulation and an estuary simulation. For the oceanic domain, the open boundary is forced with the amplitude and phase of 13 tide constituents (K1, K2, L2, M2, MU2, N2, NU2, 2N2, O1, P1, Q1, S2 and T2), extracted from the Le Provost global tide dataset (Le Provost et al., 1998). In addition, the tide potential functions of the same constituents, incorporating the effect of the earth elasticity factors, are applied to the interior of the domain. The bottom friction is expressed as a quadratic formulation of the depth averaged velocities. A constant friction coefficient Cf = 2  103 has been adopted, as this value is used in similar ADCIRC tide modelling studies, and it is equivalent to a Chezy coefficient of around 65 m1/2/s (Blain and Rogers, 1998). The horizontal eddy viscosity term mT is constant and equal to 30 m2 s1. The model is run for 10 days, with an initial 2 day long hyperbolic ramp function for minimizing computational instabilities. For the estuary model, the water elevation time series calculated in the near-shore in front of the Lima river mouth are imposed at the open boundary condition, again for a 10 day long simulation with 2 days of ramp. No river flow is imposed at the upstream boundary, while the same Cf and mT values of the outer

model are maintained. Three different simulations are performed: S5, derived from the 5 m DEM grid interpolated linearly to M5; S50, derived from the 50 m DEM grid and interpolated linearly to M50; and S5R, obtained by linear interpolation of the raw bathymetry scatter set to M5. A discussion between predicted and measured elevation and velocities at the 3 ADCPs site follows. 3.3.3.1. Free surface elevation. A visual comparison between observed and computed water elevation at sites 1 and 3 is offered in Fig. 8. Fig. 8 shows that for S5 the tide amplitudes and phases are accurately reproduced at site 1. At site 3 the model reproduces the asymmetry of the tide curve, which is apparent only in the upper part of the estuary domain, where the flow has experienced a more prolonged interaction with the bathymetry in time and space than at sites 1 and 2. At site 3 however, the model underestimates the absolute tide water elevation of around 0.5 m at high water, and of 1.0 m at low water. The vertical off-set between predictions and observations at site 3, present also for the S50 and S5R results, may be related to the model inability to meet the mass conservation law when the wet and dry option is active (Blain et al., 2010). Table 4 reports the maximum differences in the water elevations at site 3 between S5, which is taken as the reference situation, S50 and S5R, at high and low tide. Site 3 is chosen because it is the farthest from the ocean boundary, where changes between the different bathymetries might be expected to have more influence on the results. The maximum discrepancy between the three tests occurs at low water, corresponding to 2.7% of the local tide range. The fact that the differences at low water are higher than at high water may be related to a higher influence of the bathymetry over a shorter water column. The differences computed for sites 1 and 2 and not shown here, are lower than for site 3, as expected. In order to spatially represent the differences between the tests with the different kriging resolution bathymetry and the raw bathymetry, the S5, S50 and S5R water levels are plotted along a longitudinal river section in Figs. 10 and 11 for two distinct tide phases. The longitudinal river section, which corresponds approximately to the river thalweg, is plotted in Fig. 9. The depth profiles of Fig. 9 reveal that the bathymetry with a 5 m resolution tends to be identical to the raw bathymetry. The

Fig. 6. Estuary mesh M5. The navigation and main river channel show a higher resolution than nearby areas.

A.P. Falcão et al. / Journal of Hydrology 497 (2013) 152–164

159

Fig. 7. Estuary mesh M50: the spatial resolution here is generally coarser than in M5, as a consequence of less restrictive NSF values.

Fig. 8. Comparison of the S5 free surface computed elevation to observation at (a) site 1 and (b) site 3.

Table 4 Maximum difference in the water elevation between the S5, S50 and S5R simulations. Reference simulations

|S5–S50| |S5–S5R|

Water elevation difference At high water (cm)

At low water (cm)

1.2 1.1

8.9 9.6

50 m resolution bathymetry is conforming to the trend of the previous 2 bathymetries, apart from local divergences as a consequence of sudden changes of the bottom profile, visible in Fig. 9a), which are not revealed in a 50 m resolution scatterset.

Figs. 10 and 11 depict the water elevation profiles calculated respectively at the time 1am (flood tide) and 4 pm (ebb tide) of 8th October, along the same longitudinal transect. In Fig. 11 a closer view of the different profiles is introduced to make them distinguishable. Fig. 10 corresponds to an instant of the early flood tide, where the tide is entering into the estuary and propagating upstream. From the measured water elevations it can be deduced that the incoming tide has already attained sites 1 and 2 but not yet site 3, while the 3 computed profiles show an identical trend with a maximum difference of less than 5 cm. The negative gradient of the computed free surface in the upstream direction demonstrates

160

A.P. Falcão et al. / Journal of Hydrology 497 (2013) 152–164

Fig. 9. Comparison of the 5 m kriging grid, 50 m kriging grid and raw bathymetries along an estuary longitudinal transect (a) between 1000 m and 3000 m distance from the river mouth; and (b) between 5000 m and 7000 m from the river mouth.

Fig. 10. Water surface elevations for the reference bathymetries computed along a longitudinal estuary transect compared with measurements at site 1, 2 and 3 for 1 am of October 8th 2006.

the model ability in reproducing the rising tide curve along the estuary main direction. The differences, evident especially for the central measurement, are related to a slight phase lag. On the other hand Fig. 11 represents an ebb phase, with a low positive free surface gradient for both computed and measured values. The relative differences of the three water profiles are of the order of the centimeter: a close-up view of is introduced in order to distinguish the vertical gap intervening between them.

3.3.3.2. Depth averaged velocities. Comparisons between measured and computed depth averaged velocities at sites 1, 2 and 3 are plotted in Figs. 12–14 respectively.

At site 1 both velocity magnitude and phase are well reproduced for all the three simulation, with a slight overestimation of the maximum flood velocity at spring tide. The measured velocity peak occurred during October 9th may be correlated to a meteorological event. The measurements at site 2, where the instrument had to be re-deployed due to an accidental interference with the mooring, reveal a strong tide asymmetry with prolonged and intense flood phases and shorter and weaker ebb phases. An analysis of the instrument raw data showed that a strong surface flow is present and bias the averaged velocity. The model is able to capture the tide phase, but the peak flood velocity is underestimated, especially when using the 50 m kriging resolution. At site 3 the temporal succession of flood and ebb phases is well respected for

A.P. Falcão et al. / Journal of Hydrology 497 (2013) 152–164

161

Fig. 11. Water surface elevations for the reference bathymetries computed along a longitudinal estuary transect compared with measurements at site 1, 2 and 3 for 4 pm of October 8th 2006.

Fig. 12. Time series of the measured and computed depth averaged velocities at site 1. Positive values correspond to the flood phase.

Fig. 13. Time series of the measured and computed depth averaged velocities at site 2. Positive values correspond to the flood phase.

Fig. 14. Time series of the measured and computed depth averaged velocities at site 3. Positive values correspond to the flood phase.

S5 and S5, while S5R has a limited phase lag. The measured flood and ebb peaks are underestimated, with a heavy damping for S50.

3.3.3.3. Model performance indices. In order to assess quantitatively the model sensitivity to the bathymetry and mesh resolution, a

162

A.P. Falcão et al. / Journal of Hydrology 497 (2013) 152–164

Table 5 Performance indices for the three different runs, at sites 1, 2 and 3. Tide elevation MAE

Depth averaged velocity

RMSE

RMAE

NSE

MAE

RMSE

RMAE

Site 1 S5 0.095 S50 0.094 S5R 0.092

NSE

0.108 0.108 0.105

0.053 0.053 0.051

0.989 0.989 0.989

0.100 0.101 0.094

0.152 0.153 0.147

0.640 0.607 0.654

0.746 0.743 0.761

Site 2 S5 0.135 S50 0.134 S5R 0.128

0.141 0.141 0.135

0.067 0.066 0.063

0.962 0.962 0.965

0.573 0.677 0.639

0.706 0.745 0.741

0.807 0.955 0.871

0.960 1.185 1.162

Site 3 S5 0.219 S50 0.218 S5R 0.201

0.226 0.225 0.211

0.111 0.110 0.098

0.951 0.952 0.958

0.192 0.325 0.191

0.221 0.367 0.228

0.598 0.882 0.629

0.753 0.319 0.737

h P i0:5 2 obs  ycomp j; RMSE ¼ N1 N  ycomp Þ . i¼1 ðyi i i   PN obs 2 PN obs comp 2 comp obs   P yi yi  Þ  ðy y ðy yi Þ i¼1 i  = average of PN obsi¼1 2i where y RMAE ¼ N1 N i¼1  ; NSE ¼ yobs ðy  yÞ i i¼1 i observations. MAE ¼ N1

PN

obs i¼1 jyi

series of model performance indices is evaluated (Sutherland et al., 2004; French, 2010). These include the mean average error (MAE),

the root mean square error (RMSE), the relative mean average error (RMAE) and the Nash-Sutcliff coefficient (NSE). Table 5 reports the definition of these indexes and their evaluation at sites 1, 2 and 3. The range of variability of all analysed performance indices is higher for the velocities than for the water elevation, demonstrating that the model is more sensitive to the current field than to the free surface. However, the model performance is better in reproducing the water elevations: in particular excellent values are obtained at site 1 for all three model configurations, with a constant NSE equal to 0,989. In general the S5R values are better than S5 and S50. At site 2 and 3, although the absolute values are good, it is noticed a general worsening of all indices with respect to site 1: this has to be related to the vertical off-set between predictions and observations inside the estuary (Fig. 8b) and to the difficulty in reproducing the tide asymmetry, which increases with the distance from the river mouth. This trend is confirmed also for the velocity statistics, with the model performing better at site 1 than at site 3. The values of the velocity indices at site 2 are considered not indicative because of the presence of a local superficial flow which a depth averaged model is not able to reproduce. S5R shows the best MAE, RMSE and NSE at site 1, while S50 performs better for RMAE. At site 3 the best behaviour is shown by S5, while the S5R indices are, apart from MAE, worse because of a velocity phase delay (Fig. 14).

Fig. 15. Relative position of the wetted front in the central part of the Lima estuary, calculated at 4 pm of October 8th for S5 (darker line) and S5R (brighter line).

Fig. 16. Relative position of the wetted front calculated along the course of the Lima river, in the upper part of the estuary domain, at 4 pm of October 8th for S5 (darker line) and S50 (brighter line).

A.P. Falcão et al. / Journal of Hydrology 497 (2013) 152–164

3.3.3.4. Flooded area extension. Previous results are complemented by a spatial visualisation of the extension of the computed flooded area inside the estuary (Figs. 15 and 16). Fig. 15 makes a comparison between the location of the wetted fronts of the S5 and S5R simulations, calculated for 4 pm of October 8th. The wetted front indicated the margins of the inundated area of the estuary at a certain time step during the tide cycle, as a consequence of the flooding and drying of the mesh elements. The only change existing between the simulations is the use of a different bathymetry, while the mesh is the same. The maximum differences between the positions of the wetted fronts are located in the upper part of Fig. 15. These values are mainly dependent on the local bathymetry gradient, as the simulation results evidenced until show a substantial agreement. The same conclusion is valid for Fig. 16, which represents the upper part of the domain: here the higher bathymetry gradient of the river banks makes the front lines closer than in the middle estuary.

163

3.3.4. Vertical geo-referencing error In order to understand how possible errors in the vertical referencing of the bathymetry affect the computed water levels and inundated areas, the original values of the Estuary scatterset were arbitrarily increased by +0.25 m. This value is in accordance with the reported vertical offset that different field surveys, performed at different times and with different instrumentations, may experience. The new simulation, called S25+, keeps the ocean boundary conditions, model input settings, and the unstructured mesh M5 used for S5R. The results show that, for the reference station 1, no differences (less than 1 cm) are noticed, for both high and low waters, suggesting that the water surface during the dewatering process is not influenced by the new upstream bottom level. At reference station 2 the high tide water levels remain unchanged, while at low tide a maximum difference of 4.4 cm is seen, at spring tide, between S25+ and S5R. In absolute terms, the biggest differences are experienced at the most upstream station, where at low spring tide the maximum water elevation between S25+ and S5R differs of 12.1 cm, with a phase delay of around 30 min. As a consequence, the asymmetry of the tidal curve is more pronounced, which is in accordance with the enhanced interaction with the bottom morphology due to a shallower environment. No changes are instead found at high waters. The extension of the flooded areas, for the magnitude of the reported differences, do not suffer any changes; this is especially important for high water conditions, which are not affected by ordinary errors made in the vertical acquisition of the bathymetry and topography.

was considered in order to get a similar resolution as that of the original ‘‘Estuary’’ bathymetric dataset. The profile built along the thalweg channel, where field data are known in greater detail, and constructed from the 5 m cell size grid is practically coincident with the profile obtained from the raw data, obtained by merging the Estuary, Lima and contour points. Comparing with the profile built from the 50 m cell size grid, some local differences were noticed, especially in correspondence of rapidly varying slopes, as expected. However, in the simulations results, the maximum water elevation differences between the use of the 5 m and 50 m DEM grids, interpolated to the M5 and M50 finite element meshes, were of 1.2 cm at high tide and 8.9 cm at low tide. Almost the same values were obtained when the raw data was used in the modelling. The reduced model sensitivity to the use of different DEMs may be explained by the good performance of the chosen interpolator, as already said, in terms of low mean error (0.01 m). The absolute vertical elevation differences reported in this article may seem small; however these have to be related to the local bathymetric gradient and projected horizontally: for the same vertical difference, flat surfaces are inundated wider than steep slopes. This is evident in Figs. 12 and 13: where the bathymetry is flatter, the distances between wetted fronts increase, and it is on those locations where uncertainty appears. Yet, the uncertainties calculated for the terrain, not known a priori, contribute, together with other variables (quality of the model input forcings, such as river flow, tide, meteorological variables and model parameters) to make the total uncertainty in inundation studies. Regarding the influence of grid size in the flooding intertidal area, it was possible to observe that the shape of two intertidal flood areas in both simulations is consistent. In what concerns the results of the simulations when a systematic vertical difference is introduced, it was verified that the larger differences were located in upstream stations during the low spring tide, while during the high tide no significant differences were found. It can be inferred that ordinary vertical errors (0.25 m) related to possible wrong georeferencing procedures have very limited impact on the computation of the high water levels and in the flood mapping in estuarine areas. The model incapacity to reproduce the absolute values of the water levels measured in the upper part of the estuary (site 3) has to be investigated with further modelling. In particular, a test without the wet and dry option may indicate the influence of the lack of mass conservation introduced by the algorithm.

4. Conclusions

Acknowledgements

The merging of topographic and bathymetric data is frequently a requirement in coastal engineering studies, as in the case of the assessment of floods for estuarine areas. Usually the intertidal areas are not covered by the mentioned datasets, leaving a gap that prevents an accurate terrain representation, critical for hydrodynamic models and simulations. The methodology presented in this study helps to solve this problem in an accurate and inexpensive manner. A thorough selection of an interpolator and its parameters is important for the construction of an elevation surface from the input height points. Although all tested interpolators gave good validation results in terms of mean average error, ordinary kriging showed the lowest standard deviation error. The computational effort of ordinary kriging is not significantly different from that of other interpolators. Being a probabilistic estimator and considering the anisotropic behaviour of data, it was our selection of interpolator after testing several alternatives. The 5 m resolution

This work was funded by the Portuguese Foundation for the Science and the Technology (FCT), through the Doctoral Grant of Andrea Mazzolari (SFRH/BD/61161/2009). We thank the ADCIRC Coastal Circulation and Storm Surge Model Development group for providing the parallel version of the software. We want to express our gratitude to the Portuguese Army Geographical Institute, Hidrographic Institute, and the surveying company Hidrodata. We express our acknowledgment to Nuno Silva for his valuable support in fieldwork. We are also very grateful to the anonymous reviewers for their thorough reviews and constructive comments that helped to improve the manuscript. References Aguilar, F.J., Aguilar, M.A., Agüera, F., Sánchez, J., 2006. The accuracy of grid digital elevation models linearly constructed from scattered sample data. International Journal of Geographical Information Science 20 (2), 169–192.

164

A.P. Falcão et al. / Journal of Hydrology 497 (2013) 152–164

Ali, A., Zhang, H., Lemckert, C.J., 2009. Numerical of hydrodynamics of a very shallow estuarine systems – Coombabah Lake, Gold Coast, Australia. Journal of Coastal Research, Special Issue 56, 922–926. Araújo, M.A.V.C., Mazzolari, A., Trigo-Teixeira, A., 2011. Wave set-up in the modeling of storm surge at Viana do Castelo (Portugal). Journal of Coastal Research 64, 971–975. Blain, C.A., Rogers, W.E., 1998. Coastal tide prediction using the ADCIRC-2DDI hydrodynamic finite element model: model validation and sensitivity analyses in the southern North Sea/English Channel. Technical report NRL/FR/7322-989682, Naval Research Laboratory Oceanography Division, Stennis Space Center, MS, 96p. Blain, C.A., Westerink, J.J., Luettich Jr., R.A., 1994. The influence of domain size on the response characteristics of a hurricane storm surge model. Journal of Geophysical Research 99, 18467–18479. Blain, C.A., Linzell, R.S., Chu, P., Massey, C., 2010. Validation test report for the Advanced CIRCulation Model (ADCIRC) v45.11. NRL Memorandum Rep NRL/ MR/7320-10-9205 Naval Research Laboratory, Washington, DC. Buttner, O., 2007. The influence of topographic and mesh resolution in 2D hydrodynamic modelling for floodplains and urban areas. Geophysical Research Abstracts 9, 08232. Cea, L., French, J.R., 2012. Bathymetric error estimation for the calibration and validation of estuarine hydrodynamic models. Estuarine, Coastal and Shelf Science 100, 124–132. Chaplot, V., Darboux, F., Bourreannane, H., Leguedois, S., Silveira, N., Phachomphon, K., 2006. Accuracy of interpolation techniques for the derivation of digital elevation models in relation to landform types and data density. Geomorphology 77, 126–141. Cobby, D.M., Mason, D.C., Horritt, M.S., Bates, P.D., 2003. Two-dimensional hydraulic flood modelling using a finite-element mesh decomposed according to vegetation and topographic features derived from airborne scanning laser altimetry. Hydrological Processes 17, 1979–2000. Cook, A., Merwade, V., 2009. Effect of topographic data, geometric configuration and modelling approach on flood inundation mapping. Journal of Hydrology 377, 131–142. Coveney, S., Stewart Fotheringham, A., 2011. The impact of DEM data source on prediction of flooding and erosion risk due to sea-level rise. International Journal of Geographical Information Science 25 (7), 1191–1211. Dietrich, J.C., Zijlema, M., Westerink, J.J., Holthuijsen, L.H., Dawson, C., Luettich Jr., R.A., Jensen, R.E., Smith, J.M., Stelling, G.S., Stone, G.W., 2011. Modeling hurricane waves and storm surge using integrally-coupled, scalable computations. Coastal Engineering 58 (1), 45–65. DiGruttolo, N., Mohamed, A., 2011. Emerging unmanned aerial remote sending system for intertidal zone modelling: a low-cost method of collecting remote sensing data for modelling short-term effects of sea level rise, Part II: Closerange airborne remote sensing. Surveying and Land Information Science 70 (3), 119–129. Falcão, A.P., Matos, J., Casaca, J., Sousa, A.J., Gonçalves, A., 2008. Preliminary results of spatial modelling of GPS/levelling heights: a local quasi-geoid/geoid for the Lisbon area. In: Proceedings of International Symposium on Gravity, Geoid and Earth Observation, Chania, Crete, Greece, 2008. Feyen, J., Hess, K., Spargo, E., Wong, A., White, S., Sellars, J., Gill, S., 2006. Development of a continuous bathymetric/topographic unstructured coastal flooding model to study sea level rise in North Carolina. In: Spaulding M. et al. (Ed.), Proceedings of the 9th International Conference on Estuarine and Coastal Modeling, ASCE. Fortunato, A.B., Oliveira, A., Baptista, A.M., 1999. On the effect of tidal flats on the hydrodynamics of the Tagus Estuary. Oceanologica Acta 22 (1), 31–44. Fotopoulos, G., 2005. Calibration of geoid error models via a combined adjustment of ellipsoidal, orthometric and gravimetric geoid height data. Journal of Geodesy 79, 111–123. French, J.R., 2010. Critical perspectives on the evaluation and optimization of complex numerical models of estuary hydrodynamics and sediment dynamics. Earth Surface Processes and Landforms 35, 174–189. Guo, Q., Li, W., Yu, H., Alvarez, O., 2010. Effects of topographic variability and Lidar sampling density on several DEM interpolation methods. Photogrammetric Engineering and Remote Sensing 76 (6), 701–712. Hardy, R.J., Bates, P.D., Anderson, M.G., 1999. The importance of spatial resolution in hydraulic models for floodplain environments. Journal of Hydrology 216, 124–136. Horritt, M.S., Bates, P.D., Mattinson, M.J., 2006. Effects of mesh resolution and topographic representation in 2D finite volume models of shallow water fluvial flow. Journal of Hydrology 329, 306–314. Instituto Hidrográfico, n/d. Webpage of Instituto Hidrográfico (Portuguese Hydrographic Institute) . (accessed 07.12). Jones, J.E., Davies, A.M., 2007. On the sensitivity of computed higher tidal harmonics to mesh size in a finite element model. Continental Shelf Research 27, 1908– 1927. Le Provost, C., Vincent, P., 1986. Some tests of precision for a finite element model of ocean tides. Journal of Computational Physics 65, 273–291.

Le Provost, C., Lyard, F., Molines, J., Genco, M., Rabilloud, F., 1998. A hydrodynamic ocean tide model improved by assimilating a satellite altimeter derived data set. Journal of Geophysical Research 103, 5513–5529. Legleiter, C.J., Kyriakidis, P.C., 2008. Spatial prediction of river channel topography by kriging. Earth Surface Processes and Landforms 33, 841–867. Long, B., 2005. Étude hydrodynamique, sédimentologique et biologique des sites de Maria, Saint-Siméon, Bonaventure, Newport et Cap-d’Espoir dans la baie des Chaleurs, Québec, Canada. Rapport INRS-ETE pour le Ministère des Transports du Québec, 179p. Luettich, R.A., Westerink, J.J., 1999. Elemental wetting and drying in the ADCIRC hydrodynamic model: Upgrades and documentation for ADCIRC version 34.XX. Contractors Rep, US Army Corps of Engineers. Luettich, R.A., Westerink, J.J., 2004. Formulation and numerical implementation of the 2D/3D ADCIRC finite element model, version 44.XX. . (accessed 02.10.11). Luettich, R.A., Westerink, J.J., Scheffner, N.W., 1992. ADCIRC: an advanced threedimensional circulation model for shelves, coasts and estuaries, Report 1: Theory and methodology of ADCIRC-2DDI and ADCIRC-3DL. Tech Rep DRP-92-6, US Army Corps of Engineers. Mazzolari, A., Trigo-Teixeira, A., Araújo, M.A.V.C., 2012. Two dimensional unstructured mesh generation for shallow water models based on the multicriteria strategy. Actas das II Jornadas de Engenharia Hidrográfica. Lisbon, Portugal. McInnes, K., Macadam, I., Hubbert, G., O’Grady, J., 2009. A modelling approach for estimating the frequency of sea level extremes and the impact of climate change in southeast Australia. Natural Hazards 51, 115–137. Merwade, V., 2009. Effect of spatial trends on interpolation of river bathymetry. Journal of Hydrology 371, 169–181. Merwade, V., Maidment, D., Goff, J., 2006. Anisotropic considerations while interpolating river channel bathymetry. Journal of Hydrology 331, 731– 741. Merwade, V., Cook, A., Coonrod, J., 2008a. GIS techniques for creating river terrain models for hydrodynamic modeling and flood inundation mapping. Environmental Modelling and Software 23, 1300–1311. Merwade, V., Oliveira, F., Arabi, M., Edleman, S., 2008b. Uncertainty in flood inundation mapping: currents issues and future directions. Journal of Hydrologic Engineering 7, 608–620. Nahavandchi, H., Soltanpour, A., 2006. Improved determination of heights using a conversion surface by combining gravimetric quasi-geoid and GPS-levelling differences. Studia Geophysica and Geodaetica 50 (2), 165–180. Omer, C.R., Nelson, E.J., Zundel, A.K., 2003. Impact of varied data resolution on hydraulic modeling and floodplain delineation. Journal of the American Water Resources Association 39 (2), 467–475. Pe’eri, S., Long, B., 2011. LIDAR technology applied in coastal studies and management. Journal of Coastal Research, Special Issue 62, 1–5. Pinho, J.L.S., Vieira, J.M.P., 2007. Mathematical modelling of salt water intrusion in a northern Portuguese estuary. IAHS Publication 310, 277–290. Poulter, B., Halpin, P.N., 2008. Raster modelling of coastal flooding from sea-level rise. International Journal of Geographical Information Science 22 (2), 167– 182. Qiong, W., Guodong, Y., 2011. Geoid refinement of Songyuan irrigation area based on EGM2008 and GPS. International Conference Electronics, Communications and Control (ICECC 2011), 2666–2669. Raber, G.T., Jensen, J.R., Hodgson, M.E., Tullis, J.A., Davis, B.A., Berglund, J., 2007. Impact of Lidar nominal post-spacing on DEM accuracy and flood zone delineation. Photogrammetric Engineering and Remote Sensing 73 (7), 793– 804. Rebordão, I., Trigo-Teixeira, A., 2009. Tidal propagation in the Lima Estuary. Journal of Coastal Research SI 56, 1400–1404. Reguzzoni, M., Sansó, F., Venuti, G., 2005. The theory of general kriging, with applications to the determination of a local geoid. Geophysical Journal International 162 (2), 303–314. Smith, W., Sandwell, D., 1997. Global seafloor topography from satellite altimetry and ship depth soundings. Science 277, 1957–1962. Sutherland, J., Walstra, D.J.R., Chesher, T.J., van Rijn, L.C., Southgate, H.N., 2004. Evaluation of coastal area modelling systems at an estuary mouth. Coastal Engineering 51, 119–142. Wackernagel, H., 2003. Multivariate Geostatistics, third ed. Springer, Berlin. Weaver, R.J., Slinn, D.N., 2010. Influence of bathymetric fluctuations on coastal storm surge. Coastal Engineering 57, 62–70. Westerink, J.J., Luettich, R.A., Feyen, J.C., Atkinson, J.H., Dawson, C., Roberts, H.J., Powell, M.D., Dunion, J.P., Kubatko, E.J., Pourtaheri, H., 2008. A basin- to channel-scale unstructured grid hurricane storm surge model applied to Southern Louisiana. Monthly Weather Review 136 (3), 833–864. Yamano, S., 2007. The use of multi-temporal satellite images to estimate intertidal reef-flat topography. Journal of Spatial Science 52, 73–79. Zaletnyik, P., Völgyesi, L., Kirchnet, I., Paláncz, B., 2007. Combination of GPS/leveling and the gravimetric geoid by using the thin plate spline interpolation technique via finite element method. Journal of Applied Geodesy 1, 233–239.