1 Introduction

Extreme sea level (ESL) driven coastal inundation presents socioeconomic challenges for Pacific Island countries (PICs). Global mean sea-level rise is accelerating [22] and is expected to exacerbate ESL frequency and magnitude in PICs, which are expected to experience higher local relative sea-level rise than the global average [21]. This expectation is reinforced by recent observational evidence attributing much of the increase in extreme sea-level frequency since 1900 to human-driven sea level rise [6]. Despite this future threat, there are limited ESL inundation maps for PICs at island or national scales. Previous studies at these scales have analysed tide gauge records to evaluate how present-day ESL event frequencies change with local relative sea level rise (SLR). Vitousek et al. [39] observed PICs with microtidal regimes experience the largest change in ESL frequency, with present-day 50-year annual recurrence interval (ARI) events occurring every ~ 5-years on average in PICs with only 0.1 m to 0.2 m SLR. Tebaldi et al. [36] also reported that by 2100, present-day 100-year events could occur once per year on average at approximately half of global coastal locations under a 1.5 °C warming scenario, with the timing of this frequency shift dependent on the duration of sustained warming. While these findings present a significant concern for flood risk management in PICs, simulating and mapping ESL inundation is required to quantify the spatiotemporal risk of exposure and impacts to social and economic elements.

Coastal inundation models are often categorized as ‘static’ and ‘dynamic’ based on how physical processes are resolved [35]. Static models (i.e.,‘bathtub’ approach) are applied in geospatial information systems (GIS) to identify topographical elevations below an ESL elevation, with hydrological connection to the seas. These computationally efficient models support ESL inundation mapping over large geographical areas. However, efficiency comes with several inherent model limitations whereby important physical processes are not accounted such as conservation of mass for flows, landscape roughness and flood protection structure effects on floodwater spreading, and storm-driven ESLs occur over an infinite duration [32]. While dynamic models can simulate these processes, high computational expense and the limited availability of high-resolution topographical data captured by Airborne Light Detection and Ranging (LIDAR) techniques, can limit their practical application to subnational scales for many countries, including PICs. Model selection should then consider several factors, including the available water level data records (e.g., tide gauge location and duration), topographical data resolution, computational resources, and whether increasing model complexities achieve the accuracy requirements for maps to inform national-scale flood risk management decisions [20].

We present a new dataset of nationwide ESL inundation maps for Vanuatu. Our research objectives were to develop a modular ESL inundation mapping framework based on a static model approach, and to then apply the model to create inundation maps for the inhabited islands. A static model approach was chosen for several practical reasons: the limited current availability of inundation maps, limited computational resources, incomplete nationwide high-resolution topographical data, and the need to generate island- and national-scale social and economic exposure and impact assessments from ESL inundation under present-day and future higher sea levels. The model requires three primary input datasets: total water levels, island boundaries, and a digital elevation model. We describe the model framework, input data, and geoprocessing steps in Sect.  2, with ESL inundation maps interpreted in Sect.  3. In Sect.  4, we offer insights into the dataset’s limitations and its anticipated use for coastal flood risk analysis.

2 Methods

The framework presented in Fig. 1. follows a static-inclined model approach for generating inundation maps using total water level and topographic data [35]. Static-inclined models differ from a traditional static model (i.e., planar water surface along shorelines) whereby a spatially varying water surface is simulated along shorelines. Model input data represents the total water level (Sect. 2.2), island shape (Sect. 2.3) and a digital elevation model (Sect. 2.4). Several geoprocessing steps supported the input data generation and are described in the corresponding sections.

Fig. 1
Fig. 1
Full size image

ESL inundation modeling framework used in this study

2.1 Total water levels

Extreme total water level (TWL) represents the cumulative effect of several oceanographic and atmospheric processes, including mean sea level (MSL), tidal elevation, and dynamic contributions from waves and storm surge. TWL can be expressed as:

$$\:\text{T}\text{W}\text{L}\left(\text{t}\right)={\stackrel{\sim}{\eta}}_{2}+\text{I}\text{B}+\text{S}\text{L}\text{A}+\text{A}\text{T}$$
(1)

where \(\:{\stackrel{\sim}{\eta}}_{2}\) is the 2% exceedance wave induced water level nearshore based on formulation for reef-fronted coastlines [24], IB is the inverted barometer effect, SLA is the sea level anomaly, and AT is the astronomical tide. Wave runup can be a dominant component of nearshore total water level during large wave events on reef-fronted coasts, accounting for more than half of the observed water level in some cases [40], and its inclusion in the TWL definition is therefore important for capturing the primary flood driver in reef-fronted island settings such as Vanuatu. \(\:{\stackrel{\sim}{\eta}}_{2}\) was determined based on the significant wave height (\(\:Hs\)), energy wave period (\(\:{Tm}_{01}\)), and mean wave direction (\(\:{\theta\:}_{m}\)) from the ERA5 global reanalysis [17]. ERA5 has a spatial resolution of approximately 0.25° (~ 31 km) and hourly temporal resolution. IB was calculated as 1 cm of sea level rise for every 1 hPa atmospheric pressure deficit from ERA5, following the classical static inverted barometer approximation. It is noted however, this approximation can underestimate meteorologically driven sea-level responses in shallow, semi-enclosed coastal environments [26]. Sea level anomaly (SLA) was obtained from ORAS5 reanalysis [42]. For the astronomical tide, the mean high-water spring (MHWS) at each location derived from the primary tidal constituents (K1, O1, M2, S2) from the TPXO9 dataset [10], was held constant in time. These four constituents typically dominate the tidal signal across the Pacific and were considered sufficient for estimating MHWS at the national scale of this study. The additional minor constituents available in TPXO9 were not included given their comparatively small contribution to MHWS in this region. Validating the derived TWL return periods against available tide gauge records, and testing sensitivity to the full TPXO9 constituent set, represent useful directions for future refinement of this dataset. This approach helped to ensure that all extreme water level events were captured, irrespective of their coincidence with high or low tides, thereby producing conservative inundation estimates.

Given the relatively low wind setup values observed in island settings—even during tropical cyclones— the contribution of the wind to the TWL was not considered in the analysis. However, this simplification may lead to underestimation at certain locations in Vanuatu, where gently sloping island profiles and the presence of semi-enclosed embayments could significantly amplify wind setup. More broadly, the combined omission of wind setup and reliance on ERA5-derived wave and pressure forcing may underestimate the most extreme tropical cyclone-driven water levels, particularly within sheltered bays and lagoons where local wind-wave interactions are not resolved by the global reanalysis. Global storm-surge datasets that explicitly simulate wind-driven surge, such as the ERA5-forced Global Tide and Surge Model (GTSM; [25]), could provide an independent benchmark for quantifying this contribution in future iterations of this dataset, though such a comparison was beyond the scope of the present national inundation mapping approach.

TWLs were calculated at spatial intervals of 200 m to 500 m along the island coastline, depending on local geomorphology, to accurately capture variations in shoreline orientation. Waves were propagated from the closest ERA5 grid node to each coastal point with Snell’s law, which describes the bending of wave crests as they move from deep to shallow water in proportion to the local change in wave celerity [37]. TWL components were interpolated to match the 3-hour resolution of the ERA5 dataset. This interval may not resolve instantaneous wave height peaks during rapidly evolving storms, but was considered adequate for characterising the sustained wave conditions contributing to the ARI return periods analysed.

The annual maxima from the 44-year TWL time series were statistically modelled using the Generalized Extreme Value (GEV) distribution, which adapts its shape parameter to reflect regional differences in wave climate. Vanuatu experiences tropical cyclones on an annual basis [8], providing a sufficient number of extreme events within the 44-year record to support a robust annual-maxima GEV fit at each coastal point. GEV parameters were estimated at each site using maximum likelihood estimation, with confidence bounds derived from the asymptotic variance of the fitted parameters. Return level uncertainty is consequently larger at cyclone-dominated points with a positive (Frechet-type) shape parameter, where a small number of large events can strongly influence the upper tail of the fitted distribution. An annual-maxima GEV approach was adopted in preference to a peaks-over-threshold Generalized Pareto Distribution (GPD) to preserve independence between selected extremes and avoid the need to define a threshold consistently across many hundreds of coastal points nationally. While a GPD approach could make more efficient use of the available record, applying it uniformly at national scale without site-specific threshold sensitivity testing introduces its own source of uncertainty. This is consistent with experience from comparable regional studies, where GPD approaches have required extensive per-point threshold tuning to yield stable results when applied in an automated way. In regions dominated by extratropical storm swell, the GEV tends toward the Weibull distribution (bounded tail), whereas in cyclone-dominated areas, the GEV aligns more closely with Frechet or Gumbel distributions, characterized by heavier tails. By fitting the GEV at each coastal point, TWL annual return intervals (ARIs) were calculated for 5-, 10-, 25-, 50-, 100-, 250-year events.

Future rates and timing of global mean sea-level rise (SLR) remain highly uncertain, particularly beyond 2050 [13]. A scenario-based approach was adopted to account for SLR uncertainty, evaluating incremental increases of 0.1 m above present-day MSL [31]. This approach allowed the generation of a wide range of inundation scenarios independent of future SLR projections. A maximum 2 m SLR was considered to represent a range of future SLR scenarios along Pacific Island coastlines over the next century, under various Shared Socioeconomic Pathways (SSPs) [13].

TWL scenarios for each island coastline provided a proxy input for the inundation model, facilitating intersection with high-resolution digital elevation models. The methodology was designed to capture TWL gradients along the shoreline of each island, enabling the generation of two-dimensional (2D) water surface elevation maps. These spatially variable surfaces introduce a degree of complexity and realism beyond a traditional (planar), uniform static inundation mapping approach and lead to more accurate representations of potential inundation extents.

2.2 Island shape digitization

Island boundaries were extracted from the global database of administrative boundaries, GADM [16], or OpenStreetMap [28] for locations where GADM omitted islands. Extreme total water levels were calculated along the island coastline (if not reef fronted) to capture shoreline orientation and estimate the breaking significant wave height. TWL points were created using the ‘Points along geometry’ tool in QGIS, whereby the atoll rim or island ‘multipolygon’ layer was split into evenly spaced points every 200 m to 500 m. This tool also produced an ‘angle’ attribute to calculate shoreline orientation.

2.3 Digital elevation model

ESL inundation mapping requires high-resolution topographical information represented as a digital elevation model (DEM) [15]. Available DEMs in PICs include satellite-derived 30 m global DEMs [12] and 1–2 m LIDAR-derived DEMs. Two LIDAR datasets were available for Vanuatu: a 2012–2013 survey covering Efate, Santo and Malekula [14] and a 2023 survey covering Efate, Malekula and Tanna [29]. While LIDAR-derived DEMs offer higher resolution, Vanuatu has partial coverage at island or sub-island levels. A national DEM was achieved by generating consistent 5 m DEMs for inhabited islands based on down sampling LIDAR DEMs where available and up sampling the satellite-derived Forest And Building removed DEM (FABDEM V1-2) [23] where LIDAR was not available. These topographical information sources formed LIDAR, satellite and composite DEMs at island scales.

Vertical accuracy differs between the two DEM sources used in this study. LIDAR-derived DEMs typically achieve vertical accuracies on the order of 0.1–0.2 m RMSE, whereas FABDEM vertical accuracy was in the range of approximately 1–5 m RMSE, varying with land cover, terrain slope and region. Provinces without LIDAR coverage in this study (i.e., Penama and Torba) inherit greater DEM-driven uncertainty in mapped inundation extent and water depth than provinces with composite LIDAR and FABDEM coverage, particularly across low-lying coastal land where small elevation errors can translate into large errors in inundation extent and water depth.

Vertical datum consistency between the DEM and the TWL outputs used to derive inundation surfaces is a further source of uncertainty [15]. LIDAR and FABDEM elevations used in this study are referenced to MSL through EGM2008 geoid for the FABDEM and 2020-era LiDAR of Efate, Malekula and Tanna, and through a local tidal datum framework for the 2010-era LiDAR of Efate, Santo and Malekula, while TWL are referenced to local mean sea level (MSL) as represented in the ORAS5, TPXO9 and ERA5 products used to derive its components. The offset between the geoid and the local time-mean sea surface (mean dynamic topography) can reach several tens of centimetres in parts of the Pacific [34], consistent with cross-island vertical datum offsets of similar magnitude documented elsewhere in the western Pacific (e.g., 42–54 cm between Taiwan and its offshore islands; [19]) Any unresolved offset between the different MSL approximating vertical references represents an additional source of uncertainty in the mapped inundation extents and water depths, alongside the DEM vertical accuracy uncertainty discussed above.

2.4 Inundation mapping

ESL inundation was mapped using a static-inclined approach as described by Stephens et al. [35]. This has several advantages for nationwide inundation mapping, which include: (1) implementation in a GIS environment using geoprocessing and spatial interpolation functions; (2) low computational demand for operation on a standard personal computer. Inundation was calculated for each island using the Generic Mapping Tools (GMT) to interpolate an inland water surface from TWL points at the resolution of the 5 m DEM. From the TWL points, a nearest-neighbour interpolation was first used to fill the inland water surface at 100 m resolution (GMT `grdfill`), which was then downscaled to the 5 m DEM resolution using bilinear interpolation (GMT `grdsample`). Inundation grid cells were retained when at least one of their cardinal neighbors were inundated and were hydraulically connected to the coastline [41]. Water depth above ground was computed for DEM grid cells from the difference between the inland water surface and the underlying terrain elevation.

2.5 Validation

The ESL static inundation model approach was validated against a hydrodynamic model for Luganville township, Santo Province, Vanuatu [1]. The hydrodynamic model was developed independently from the static-inclined model presented in this study, but enables for spatial comparisons of inundation extent and water depth between the two datasets. The hydrodynamic model used the GPU-optimised, wave-hydrodynamic coupled model XBGPU [3], a variant of XBeach [33], simulated at 10 m resolution for the Luganville township domain. Boundary conditions were derived from the Van-KIRAP probabilistic tropical cyclone modelling system [18, 38], which combines a 10,000-year synthetic tropical cyclone track catalogue with a joint-probability analysis of offshore significant wave height and storm-tide level for each ARI. The wave-driven component of the boundary water level was derived using the same Merrifield et al. [24] reef-fronted formulation adopted for η2 in this study (Sect. 2.2). The model was calibrated and validated against the Luganville tide gauge record during Tropical Cyclone Harold (April 2020).

A spatial comparison for ARI100 ESL inundation is presented in Fig. 2, showing: (a) the static model derived ARI100 ESL inundation extent; (b) the hydrodynamic model derived ARI100 ESL inundation extent; (c) the symmetric difference between modelled ARI100 ESL inundation extents; and (d) the spatial distribution of mean absolute error (MAE) in modelled ARI100 ESL water depth. In Luganville, the static model produced an ARI100 ESL inundation extent of 47.8 ha compared to 20.9 ha for the hydrodynamic model. The two models shared an overlapping inundation area of 20.7 ha (43.4%). The symmetric difference totalled 27.4 ha, dominated by static-only inundation (27.1 ha; 99%) with minimal hydrodynamic-only area (0.3 ha; 1%). This asymmetry reflects the tendency of static inundation model approaches to overestimate inundation extent in the absence of flow resistance constraints [9, 32]. Within the overlapping inundation area, water depth MAE was 0.83 m (RMSE: 1.42 m) with a bias of − 0.45 m confirming the hydrodynamic model predicts lower water depths than the static model. Higher water depth MAE values are observed along river channels and low-lying coastal plains where the hydrodynamic model resolves flow routing and bathymetric effects. Lower water depth MAE values were observed along steeper, protected shorelines.

Some discrepancies between the static and hydrodynamic inundation extents and water depths may arise from differences in the boundary extreme sea-level inputs used by each model, rather than from the inundation modelling approach alone. Both models draw on the same underlying Merrifield et al. [24] and Becker et al. (2014) reef-fronted formulation for the wave-driven component of TWL, so the inclusion of this component is not itself a point of methodological difference between the two model approaches. However, the two models derive TWL from different data sources and extreme-value methodologies. This study fits a GEV distribution to annual maxima from a 44-year ERA5 reanalysis-based TWL hindcast applied nationally, whereas Andrews et al. [1] derive boundary conditions from a 10,000-year synthetic tropical cyclone track catalogue (Van-KIRAP; [18]) with joint-probability analysis of offshore wave height and storm-tide level, specific to southeast Espiritu Santo. Notably, Andrews et al. [1] found that their probabilistic synthetic-cyclone wave heights were substantially higher at equivalent ARIs than an earlier 40-year-hindcast-based study for the same location [4]. This hindcast approach is methodologically similar to the one applied in this study, suggesting the boundary TWL inputs for the two inundation models compared in this validation may differ independently of the modelling approach itself. A further, likely more fundamental, contributor is that the two modelling approaches treat energy dissipation differently: the static-inclined model floods all land below the total water level without attenuation, whereas XBGPU explicitly resolves the reduction in flow energy and water level as waves and currents propagate across the reef flat, beach and inland terrain. This dissipation offers a plausible explanation for why the hydrodynamic model produces less inundation despite Andrews et al. [1] reporting higher boundary wave heights than an equivalent hindcast-based approach would suggest, and reinforces reef and terrestrial dissipation as a first-order process to resolve in any future comparison using a common extreme sea-level input dataset. Testing and comparison of static and hydrodynamic modelling approaches in Vanuatu using a common, comparable extreme sea-level input dataset is needed to isolate the contribution of model physics (e.g., flow resistance, reef dissipation) from that of boundary forcing to the observed differences in inundation extent and water depth, and represents an important direction for future work. A distance-based attenuation scheme (e.g., Enriquez et al. [11] could also be explored as a complementary refinement to reduce static-model overestimation without requiring full hydrodynamic modelling nationally.

3 Data interpretation and intended use

The Vanuatu national extreme sea-level (ESL) inundation maps provide spatially explicit representations of potential coastal flooding under present-day and future sea-level rise scenarios. Each map is generated at a 5 m grid resolution, ensuring sufficient detail for island-scale exposure analysis while maintaining computational efficiency for nationwide coverage. All maps are projected in the Universal Transverse Mercator (UTM) Zone 58 South coordinate system (EPSG:32758), which is the standard for Vanuatu, enabling seamless integration with other national geospatial datasets.

The inundation maps depict areas of land inundation and water depth above ground for six annual recurrence intervals (ARI: 5, 10, 25, 50, 100, and 250 years) and incremental sea-level rise scenarios of 0.1 m, between 0 m and 2 m above present mean sea level. These scenarios provide flexibility for assessing both short-term and long-term coastal hazard exposure under different climate change pathways.

Fig. 2
Fig. 2
Full size image

A spatial comparison of ARI 100 ESL inundation extent and water depth between the static and hydrodynamic models at Luganville, Santo Province, Vanuatu. The four-panel figure shows: (a) the static model derived inundation extent; (b) the hydrodynamic model derived inundation extent; (c) the symmetric difference between modelled inundation extents; and (d) the spatial distribution of mean absolute error (MAE) in modelled water depth

Maps for provinces with LIDAR coverage (Shefa, Sanma, Malampa, and Tafea) offer higher confidence in inundation extent and depth due to the use of composite DEMs combining LIDAR and FABDEM data. Provinces without LIDAR coverage (Penama and Torba) rely on satellite-derived DEMs, which introduce greater uncertainty, particularly in low-lying coastal areas where small elevation errors can significantly affect inundation predictions.

All inundation maps are available in GeoTIFF (.tif) format via the Zenodo platform [30], supporting integration into GIS workflows for exposure and risk assessments. These datasets can enable national and provincial agencies to identify areas most vulnerable to ESL-driven inundation, in order to prioritise local high-resolution inundation and risk modelling suitable to inform local climate adaptation initiatives.

Despite the static-inclined model tendency to overestimate inundation extent and water depth relative to a hydrodynamic model benchmark (Sect.  2.5), this conservative bias supports rather than undermines the dataset’s primary intended use in national and subnational screening of coastal flood risk. Consistent over-representation rather than under-representation of inundation areas, means the dataset provides a first-order basis for identifying coastal flood risk hotspots. Across Vanuatu, areas susceptible to inundation from multiple ARI and lower SLR scenarios could be identified as a higher priority for higher-resolution hydrodynamic modelling to inform adaptation interventions such as coastal protection structures, managed retreat planning, or building design standards, whereas areas exposed to fewer ARIs and higher SLR scenarios could be determined as lower priority with reasonable confidence. This screening role is the intended application of the dataset and is consistent with its use for province- and island-level hazard-exposure and risk modelling to inform adaptation prioritisation nationally.

4 Limitations and outlook

The ESL inundation maps presented in this study provide a first-order dataset that can support national- or island-level coastal hazard exposure and risk analysis under present-day and future sea-level rise scenarios. The maps extend on existing local inundation maps [7] for specific ARI scenarios (e.g., 100-year), providing a broader event-set of ARI scenarios for nationwide coastal flood risk analysis. The maps do not, however, capture fine-scale local hydrodynamic and geomorphic processes often needed for detailed, local-scale, engineering applications. Global wind and wave models perform well in describing average wind and wave climate but are known to dramatically underpredict extreme wind and wave events such as tropical cyclones (e.g., [2, 27]). Extreme cyclone conditions at higher annual recurrence intervals are therefore likely underestimated in this study. The static inundation modeling approach offers computational efficiency but simplifies the hydrodynamic processes driving ESL inundation in Pacific Island countries when compared to dynamic models [5]. This simplification is likely to overestimate inundation extents and depths, particularly for Vanuatu’s northern atolls and low-lying land and where LIDAR DEMs were not available. Nevertheless, the modular inundation mapping framework facilitates iterative map improvements at the island level, whereby low-resolution, satellite-derived DEMs can be replaced with high-resolution LIDAR, and total water levels (TWL) can be updated periodically or after major extreme sea level events (e.g., tropical cyclones). In addition, as dynamic inundation modelling methods are developed and applied in Vanuatu, the derived maps can replace static inundation maps, especially for locations with high social and economic value.