A modest 0.5-m rise in sea level will double the tsunami hazard in Macau

See allHide authors and affiliations

Science Advances  15 Aug 2018:
Vol. 4, no. 8, eaat1180
DOI: 10.1126/sciadv.aat1180


Rising sea levels will have overwhelmingly negative impacts on coastal communities globally. With previous research focused on how sea-level rise (SLR) affects storm-induced flooding, we show that SLR will also increase both the frequency and the intensity of tsunami-induced flooding, another significant coastal hazard associated with sea-level extremes. We developed probabilistic tsunami inundation maps for Macau, a densely populated coastal city located in the South China Sea, under current sea-level, 0.5-m SLR, and 1-m SLR conditions, using an extensive Monte Carlo tsunami inundation simulation. Our results indicate that conservative amounts of SLR of 0.5 m (by 2060) and 1 m (by 2100) would dramatically increase the frequency of tsunami-induced flooding incidences by a factor of 1.2 to 2.4 and 1.5 to 4.7, respectively.


Intensified coastal flooding hazard is a key consequence of sea-level rise (SLR) under future scenarios of global climate change (1). Previous studies have shown that SLR is expected to become the dominant driver in increasing coastal flooding during storm surges (13). As sea level rises, there is a resultant increase in coastal flooding hazard coincident with more frequent storm surge–induced overtopping by relatively smaller events (4).

The Macau Special Administrative Region, China, is a low-lying, coastal city state located on the northern coast of the South China Sea, a region vulnerable to SLR and extreme coastal flooding events from storms (5) and tsunamis (6). Storms and tsunamis present a significant hazard owing to population density (~21,400 people/km2) (7), a rapidly developing economy, and areas of low-elevation reclaimed land (>60% of Macau is reclaimed) (8, 9). While the threat of storm surge in this region is well known (5) and was graphically evidenced by Typhoon Hato in late 2017 (10), the threat of a tsunami hazard is not. The tsunami threat in the South China Sea primarily comes from megathrust earthquakes generated along the Manila Trench (Fig. 1A). There are multiple lines of evidence suggesting that the Manila megathrust is capable of generating basin-wide tsunamis (11, 12). These include infrequent large earthquakes, high convergence rate between the Philippine Sea Plate and the Sunda Plate, geological records, and historical records (fig. S1 shows the tsunami records in the South China Sea, and Materials and Methods outlines our geophysical records). Such infrequent yet potentially devastating tsunami hazard has rarely been taken into account in this region’s preparedness efforts, such as developing tsunami evacuation routes or adjusting building codes (13). While some assessments of tsunami hazard exist (6, 14), how the future tsunami hazard evolves with a rising sea level has not been studied.

Fig. 1 The South China Sea region with nested simulation domains.

Rectangular boxes outline the nested simulation domains. (A) The South China Sea and the three zones of the Manila megathrust. The convergence rate and coupling ratio are marked in white background boxes. An example of slip distribution from a Mw 8.0 earthquake in zone III is shown. (B) The geographical location of the Pearl River Estuary and major cities. (C) The elevation map of Macau overlaid by the land-reclamation history modified from (51). Note that the artificial islands A, B, and E2 are still under construction, while C, D, and E1 are planned reclamation.

To address this knowledge gap, we first prepare probabilistic tsunami inundation maps based on a 100,000-year synthetic earthquake catalog for the Manila megathrust (Table 1). The catalog is primarily estimated from a geodetic coupling model (15) by assuming that slip potency accumulated on the megathrust will be released by seismicity. Assuming that the seismicity follows the truncated Gutenberg and Richter (G-R) relationship and allowing a coupling ratio uncertainty of ±50% (12), we subdivide the Manila megathrust into three zones and follow the same approach used by Li et al. (6) to generate earthquakes with heterogeneous slip distribution (fig. S2) and randomly choose a rupture area from a predefined rupture area database (see Materials and Methods).

Table 1 Earthquake return period and earthquake number in a 100,000-year synthetic catalog.
View this table:

We use COMCOT (Cornell Multi-grid Coupled Tsunami Model) (16) to simulate the tsunami inundation process for all tsunamigenic earthquake events with a moment magnitude (Mw) greater than or equal to 8.0 in the synthetic catalog (see Materials and Methods). In total, we simulate tsunamis of 1886 synthetic earthquakes. We use four nested grids (Fig. 1), with the bathymetric and topographic grid resolution varying from ~2000 m in the source region to 25 m in the city area. Nautical charts with scales ranging from 1:5000 to 1:250,000 in the Pearl River Estuary near Macau and 4-m resolution topographic data (measured by the Macau Cartography and Cadastre Bureau) are used to produce inundation maps with a high degree of detail. For each earthquake magnitude, up to 821 scenarios with diverse stochastic source models are calculated to account for the variability of complex earthquake sources (17). Simulation results from all the 1886 scenarios are then aggregated at each inland location to calculate their hazard level (see details in Materials and Methods).

To investigate the influence of SLR on the hazard assessment results, we recalculate all the synthetic tsunami events at two projected SLRs for Macau (18). On the basis of output data of multiclimate models derived from the Coupled Model Intercomparison Project Phase 5 (CMIP5) and records from tide gauge and satellite altimetry, the SLR rate in Macau is estimated to accelerate from 1.35 mm/year during 1925–2010 to 4.2 mm/year over 1970–2010 (18). These values are higher than the global mean SLR rate of 3.2 ± 0.4 mm/year for the same period (19, 20). On the basis of projected local SLRs of 30 to 51 cm by 2060 and 65 to 118 cm by 2100 (18), we choose 0.5- and 1-m SLR values to represent the sea levels by the mid-century and at end of this century, which we deem as the physically plausible scenarios (21, 22). Our choice was made on the basis of regional and global SLR studies and the evolved understanding of sea-level projections since the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (23). With the simulations of an additional 1886 events under 0.5-m conditions and another 1886 events under 1-m SLR conditions, we derive three sets of probabilistic inundation maps in total.


Probabilistic tsunami inundation scenarios at current sea level

Probabilistic tsunami inundation scenarios at current sea level (Fig. 2) indicate that the highest tsunami hazard lies in the more exposed coastal areas in southern Coloane (P5 and P6; see Fig. 2C for map locations of sites P1–P6). Notably, the densely populated Macau Peninsula, especially the inner harbor area (P1) in the west, is the most susceptible area to tsunami hazard because of its relatively low elevation (less than 2 m above mean sea level). In contrast, only limited portions of coastlines in Taipa and Cotai are hazardous even when the return intervals exceed 500 years. Although reclaimed and low-lying, the densely populated and economically important Cotai strip is protected by the rocky Taipa Island and Coloane Island in the north and south. We note that at contemporary sea level, the potential tsunami impact at both 200- and 500-year return periods is minimal with inundation limited to a narrow low-lying area of the inner harbor (Fig. 2, A and B). Higher inundation depths and flood extents at current sea level only occur at longer return periods, exceeding 1000 years (Fig. 2C). At 1000-year return intervals, the expected maximum water depth reaches 3 to 5 m in the southern Coloane coastal area (P5 and P6) and up to 2 m in the western (P1 and P2) and northeastern Macau Peninsula (P3) where the inundation distance exceeds 1 km inland.

Fig. 2 Tsunami inundation maps for different sea-level conditions.

(A to C) Current sea level, (D to F) SLR 0.5 m, and (G to I) SLR 1 m with 200-year (left column), 500-year (middle column), and 1000-year return periods (right column).

Placed in the context of earthquake hazard and assuming that the 50th percentile of inundation heights among all simulations for a given earthquake magnitude represents the average impact of the corresponding earthquake magnitude, our data suggest that, with current sea level, tsunamis produced by earthquakes with a Mw of <8.4 or even 8.6 would have very limited impact on contemporary Macau (fig. S3, A to C). Our modeling shows that perceptible inundation (>1 m) only occurs when earthquakes exceed a Mw of 8.6 and is confined to the southern Coloane area (P5 and P6; fig. S3C). It is only when earthquakes reach a Mw of 8.8 that large areas of the Macau Peninsula start receiving >1-m inundation, with corresponding inundation depths of >5 m in the southern edge of Coloane (fig. S3D). However, at the top end of possible Manila Trench earthquakes, results do show that large earthquakes (Mw of ~9.0) could cause severe inundation and considerable damage in the greater part of Macau with >2-m inundation depth, including the sheltered Cotai strip area.

While we note that the 50th percentile of the inundation heights only represents the median impact of a specific earthquake magnitude, we acknowledge that the impact can be very diverse owing to source complexity (6). This diversity was shown using a suite of tsunami scenarios produced by earthquakes with a Mw of 8.6 and 8.8 (fig. S4), where the difference between the 10th and 90th percentiles of inundation depth could vary as much as 1 to 2 m. The driver of the observed diversity is likely the randomly selected rupture locations and stochastic slip distribution model (6). Our rupture model extends from the southern offshore of Taiwan to the western offshore of Palawan, and earthquakes that occur in the northern portion of the Manila Trench have a much larger tsunami impact on Macau owing to the source directivity and bathymetry effects (24). The sensitivity of near-field tsunami inundation to earthquake rupture complexity has also been highlighted by many recent studies (17, 25), and by definition, Macau does not belong to the near-source field since the tsunami travel time from the Manila Trench to Macau is more than 3 hours. Here, our work suggests that the heterogeneous feature of slip distribution still plays a significant role in determining the tsunami impact in the relatively far field.

The effect of SLR

To examine the effect of SLR, we produced expected inundation maps in 200-, 500-, and 1000-year return periods at current sea level (baseline, Fig. 2, A to C), 0.5-m SLR (Fig. 2, D to F), and 1-m SLR (Fig. 2, G to I). With an SLR of 0.5 m, the inundation extent of the 500-year return interval quickly expands in the low-lying inner harbor area (P1; Fig. 2E); the larger area in the Macau Peninsula could be inundated with up to 2-m water depth at the 1000-year return interval (P2 and P3; Fig. 2F). If sea level is raised 1 m, then the tsunami impact with the 200-year return period (Fig. 2G) exceeds that of the 500-year return period at current sea level (Fig. 2C). The inundation map with the 500-year return period (Fig. 2G) is similar to that with the 1000-year period at current sea level (Fig. 2C). With 1 m of SLR, the 1000-year return period tsunami exceeds that of the majority of Mw 8.8 events at current sea level (fig. S3D).

We generated hazard curves at six representative locations (P1 to P6; these coastal locations are either directly exposed to the sea or sheltered) to compare the inundation depths with different sea-level conditions (Fig. 3), and we find that, with 0.5 and 1 m increased sea level, the inundation depths increase significantly, although the increased values vary across the six locations (table S1). For example, at location P1, the median inundation depth increases from 1 m at current sea level to 1.5 m at 0.5-m SLR and 1.8 m at 1-m SLR at the 1000-year return period. At location P6, the inundation depth of the 1000-year return period increases from 1.3 to 2.2 m and 2.8 m at these three sea levels (see table S1). It should be noted that the increase of inundation depth can exceed the increase of sea level. In some locations (for example, P3, P5, and P6), the inundation depths may be 0.2 to 0.5 m greater than the 0.5- and 1-m increased sea level. The underlying driver of this magnified flooding is the higher background water level (deeper water depth) from which tsunami wave can build upon. Similar SLR-induced amplification in flood potential is reported by previous studies [for example, (3)], which demonstrated that SLR results in storm tides with greater amplitudes and higher run-up due to depth and friction change. We emphasize here that the impact of SLR in the coastal region is a nonlinear process (for example, the energy dissipation due to bottom friction, the interaction between land and waves, etc.), which can only be better addressed by taking other processes into account rather than just amplifying the baseline case by an equivalent SLR.

Fig. 3 Hazard curves at six selected locations for different sea-level conditions.

Black, blue, and red solid lines show hazard curves for current sea level, 0.5-m SLR, and 1-m SLR, respectively. See locations in Fig. 2C. (A) Inner harbor P1, (B) Kai Chi Kei P2, (C) northeast Macau Peninsula P3, (D) west coast of Cotai P4 (E), southwest coast of Coloane P5, and (F) southern coast of Coloane P6. The black, blue, and red shaded areas depict the hazard curves produced by 50% higher earthquake frequency and 50% lower earthquake frequency (representing the upper and lower bounds of hazard levels) under current sea level, 0.5-m SLR, and 1-m SLR, respectively.

We note the wavy character of these hazard curves, especially around the return period of 2324 years, which corresponds to Mw 9.0 earthquakes. The reason for this observation is likely twofold: First, tsunami wave heights generated by Mw 9.0 earthquakes are much higher than those generated by Mw 8.8 earthquakes. Second, larger earthquake magnitudes exhibit much longer return periods. The wavy character could also be caused by the wide variety of uncertainties in the fault parameters. Given the limited knowledge about large subduction zone earthquakes (for example, the rupture dimension for Mw ≥8.6 events), it is very difficult to incorporate an appropriate uncertainty associated with the rupture length and width in Probabilistic Tsunami Hazard Assessment (PTHA) studies. Our data set for return periods of less than 2000 years is sufficiently large to assume that the respective results are very robust.

We calculate the frequency change of these locations experiencing the same water depth by dividing the return period at current sea level with the one at 0.5- and 1-m SLRs (fig. S5). We observe that 0.5- and 1-m SLRs would increase the frequency of tsunami-induced flooding incidences by a factor of 1.2 to 2.4 and 1.5 to 4.7, respectively. For example, in the northeast of Macau Peninsula P1, the return periods of experiencing 1-m inundation depth under current sea-level, 0.5-m SLR, and 1-m SLR conditions are 962, 800, and 458 years, respectively (Fig. 3A). The return periods of experiencing 1-m water depth at location P5 are 407, 179, and 129 years (Fig. 3E). Significant reductions in return period of 1-m tsunami flooding can also be seen in other selected locations (Fig. 3 and table S1). The significantly increased frequency is contributed mostly by earthquakes with a smaller magnitude. As shown in Fig. 4, the majority of earthquake scenarios with Mw 8.2 to 8.4 do not cause tsunami inundation at current sea level (Fig. 4, A and B). In contrast, with 1-m SLR, the same group of events produces moderate to significant flooding by increasing the inundation depth up to 1 m and expanding the inundation extent in the inner harbor area (P1; Fig. 4, D and E). As sea level rises, Manila Trench earthquakes that do not pose a tsunami threat at current sea level become potentially damaging as the same tsunami waves combined with higher sea level lead to inundation that is deeper and extends further inland. The implication is that, Macau, though relatively safe from tsunami today, will rapidly become tsunami-prone with SLR. The same is likely true of many other city states (for example, Hong Kong), densely populated island communities (for example, Malé, Maldives), or burgeoning urban archipelagos (26).

Fig. 4 Intensified tsunami inundation produced by same magnitude earthquakes under 1-m SLR conditions.

For Mw 8.2, 8.4, and 8.6, (A), (B), and (C) show the 90th percentile of the inundation depth for current sea level, and (D), (E), and (F) show the 90th percentile of the corresponding inundation depth for 1-m SLR.

Synthetic earthquake frequency

The limited observation time span of seismic data along with the spatial coverage of geodetic data means considerable uncertainties exist in the coupling ratio between the Sunda Plate and the Philippine Sea Plate. To account for this, we considered a conservative range of megathrust coupling rates along the trench from a lower-bound 0.12 to 0.15 coupling ratio to an upper-bound 0.38 to 0.45 coupling ratio. The range of coupling ratio falls in the estimation provided by Hsu et al. (12) who derive coupling rates of 0.34 or 0.48 along the Manila subduction zone at latitudes 15° to 19°N using two different block models constrained by the Global Positioning System (GPS) data collected from 1998 to 2015 in the Luzon Island.

We produced hazard curves for Macau using the two end-member scenarios, with the baseline scenario (Fig. 3) demonstrating that earthquake frequency has a significant impact on the hazard level as it shifts the hazard curve toward shorter recurrence intervals for any given inundation depth and higher (deeper) inundation for a given return period. The hazard curves also indicate that for the 0.5-m SLR condition, the curves coalesce for return periods <1500 years, suggesting that the increased tsunami hazard caused by 0.5-m SLR in Macau is more or less equivalent to a 50% increase in earthquake frequency along the Manila Trench. This similarity can also be observed from the inundation maps generated under these two conditions (compare Fig. 2, E and F, with fig. S6, G and H).

Worst-case scenario

The three-zone segmented model applied thus far does not include large rupture-across-zone (RAZ) earthquakes that could rupture through the whole Manila Trench. Several RAZ earthquakes have occurred in the past, including the 2011 Tohoku-oki earthquake rupturing through large segments of the Japan Trench (27), the 1868 Arica earthquake rupturing across the Nazca Ridge (28), and the 2007 Solomon Islands earthquake rupturing across a triple junction in the western Pacific (29). We address this issue with a worst-case scenario of a Mw 9.3 earthquake in which faults rupture through the rugged subducting seamounts (fig. S7A) that was modified from the fault model of (11), which assumed a full coupling between the Philippine Sea Plate and the Sunda/Eurasian Plate.

The inundation map for the worst-case scenario at current sea level (fig. S7B) shows that almost all reclaimed land in Macau would be flooded (fig. S7B), leaving only the original islands untouched. Water depths exceeding 5 m would be experienced in the inner harbor, the northeast Macau Peninsula, and the southern Coloane Island. The southern Macau Peninsula and the two artificial islands could expect 3- to 5-m water depth, while the water depth in the western portion of the Cotai strip and the University of Macau could reach up to 3 m.

The probability of a 1-m SLR by 2100

While we demonstrate the potential impact of 0.5- and 1-m SLRs in significantly enhancing the tsunami hazard in Macau, we fear that the projected 1-m SLR by the end of this century could be a conservative estimate for this region. The key projection data in our study (18) are obtained from CMIP5, known for insufficient consideration of projections of the Antarctic ice sheets (30). New observations including the sustained increase in ice discharge from the Amundsen Sea Embayment, West Antarctica, during the last few decades (31) and the widespread grounding line retreat of some islands in West Antarctica (32) along with advanced modeling imply that the contribution of Antarctica to future SLR could be much larger than previously thought (33). The coupled model, which links atmospheric warming with the dynamic change of the Antarctic ice sheets, suggests that Antarctica alone has the potential to contribute more than a meter of SLR by 2100 globally (33). On top of the projected SLR caused mostly by human activity, recent work at Belitung Island, Indonesia, together with observations along the south coast of China recorded up to 60 cm of sea-level fluctuation on the centurial time scale nearly 7000 years ago that underscores the potential for more than half-meter sea-level fluctuation in the region without the influence of human-induced climate change (34). It is therefore feasible that the region may experience a ~1.5-m SLR by 2100. We address this possibility with an analysis of a 1.5-m SLR scenario, as presented in the Supplementary Materials (fig. S8).

The effect of SLR on the reclaimed islands

Increasing sea level and the resultant increase in tsunami hazard raise concerns for Macau’s extensive land reclamation plans. The likelihood of newly reclaimed lands resisting extreme future flooding events under SLR depends largely on the initial elevation designed for these new lands. We quantitatively tested this idea by assuming two different initial elevations: 3 and 5 m for the artificial islands A, B, C, D, and E1 (Fig. 5). We simulated 108 Mw 8.6 earthquakes under current sea-level, 0.5-m SLR, and 1-m SLR conditions and produced inundation maps (Fig. 5) that show that with 0.5-m SLR and land heights of 3 m, islands A and B would easily be overtopped by tsunamis from Mw 8.6 earthquakes, while all new islands would be inundated with up to 2-m water depth with 1-m SLR and Mw 8.6 earthquakes. In contrast, a design elevation of 5 m would allow most islands to sustain tsunamis generated by Mw 8.6 earthquakes with a return period of around 1000 years that can also accommodate 1 m of SLR. This simple test underscores the importance of considering SLR in the design of future land reclamations and implies that similar approaches should be expanded to designs of other coastal infrastructure potentially threatened by rising sea levels.

Fig. 5 Higher elevation of newly reclaimed islands crucial for preventing future tsunami inundation.

Comparisons of the 90th percentile produced by 108 individual inundation maps based on Mw 8.6 earthquakes with 3-m initial elevation for (A) current sea level, (B) 0.5-m SLR, and (C) 1-m SLR, and 5-m initial elevation for (D) current sea level, (E) 0.5-m SLR, and (F) 1-m SLR. Note that the slightly increased inundation in the northern Macau area (D to F) is caused by wave reflection between the higher-elevation reclaimed islands and the mainland.


Using Macau as a case study, we have demonstrated that modest levels of 0.5- and 1-m SLR significantly increase tsunami hazard in Macau and shorten the return period of potential tsunami flooding for most of the territory. Notably, rising seas increase the tsunami impact generated by smaller earthquakes, which present little threat of inundation under current conditions.

Our study shows that tsunamis, like other sea-level extremes, for example, storm surges and extreme tides, will likely inundate low-lying coastal communities more frequently and more destructively as sea level rises. In contrast to storm surge events, tsunami events commonly have a much larger spatial footprint that is regional or even basin-wide if the source magnitude is large enough (for example, the 2004 Indian Ocean tsunami and the 2011 Tohoku-oki tsunami). Although the vulnerability of coastal communities to rising seas is well known, the capacity for a single tsunami event to devastate coastal areas is rarely considered in coastal planning.

We acknowledge that our study has several limitations in that we only consider the mean high water (MHW) as the tidal stage for all simulations (see details in Materials and Methods). Where the tide range is large, it may be pertinent to incorporate the influence of tides on the tsunami inundation in a probabilistic way. This would require an extension of this study and a new approach. We also simplified the overland flow by not including buildings and hence using the bare-ground data of Macau. As such, we cannot account for complex flow patterns (for example, wavefront collisions with buildings) and channeling effects, which likely cause local increases in water depth and velocity. Including the city buildings in the probabilistic tsunami inundation simulations is highly computationally expensive and well beyond the scope of this study. Nevertheless, together with coastal population and infrastructure data, the probabilistic inundation maps provided in this study can be applied to calculate tsunami exposure and risk, which can guide tsunami risk reduction efforts.

Tsunami hazard with its low recurrence but high spatial footprint is commonly overlooked in future planning; however, it is clear that for coastal cities such as Macau and others that are reliant on reclaiming coastal land, an integrated approach to tsunami hazard analysis that fully considers rising seas is essential for developing sustainable coastal infrastructure for the coming centuries.


Some of the following methods are similar to those previously published (6).

Determination of earthquake return period along the Manila Trench

Megathrust earthquakes generated along the Manila Trench are recognized as the biggest tsunami threat in the South China Sea (6, 11, 13, 14). Of particular concern are the contrasting facts that no earthquake larger than Mw 7.8 has occurred along the trench since the Spanish colonization of Luzon in the 1560s (35) and the fast convergence rate (83 to 94 mm/year) between the Sunda Plate and the Philippine Sea Plate (15), which is comparable to other well-known giant earthquake–generating megathrusts (for example, Japan, Alaska, Sunda, and northern Chilean megathrusts) (36).

Table 1 shows the a and b values of the G-R relationships of the synthetic catalog and the discrete return periods for different magnitude spans. The synthetic catalog was generated largely from the first-order GPS-based megathrust coupling model from Hsu et al. (15) and using the approach suggested by Ader et al. (37) and the truncated G-R relationship. As both the plate convergence rate across the Manila megathrust and the average megathrust coupling ratio change from north to south along the megathrust (15), we further separated the entire megathrust to three different zones (zones I, II, and III; see Fig. 1A) to better capture the convergence rate at different parts of the megathrust fault (Fig. 2A). The boundaries of these three zones are informed by the presence of the Scarborough seamount chain between latitudes 15° and 18°N, where the oceanic floor is highly fractured. We assumed that the maximum magnitudes of earthquakes in each zone were constrained by the length and width of that zone and were Mw 8.4, Mw 8.4, and Mw 9.0 for zones I, II, and III, respectively.

Among the three zones, an open question is whether the northern portion of the Manila subduction zone (zone III) has the capacity to produce Mw 9.0 earthquakes. Because of the limited spatial and temporal coverage of the observational data, large uncertainty exists in the long-term seismogenic behavior of this zone. However, several geological features indicate that the northern Manila subduction zone is likely to be able to generate large earthquakes: (i) The fault geometry obtained from earthquake focal mechanisms suggests that the dip angle in zone III ranges between 5° and 15°, which is much flatter compared to the steeper dip angle in zone I. According to Bletery et al. (38), mega-earthquakes preferentially rupture on flat megathrusts. The dip angle in zone III is in a similar range with those found in other subduction zones, for example, the Japan-Kuril-Kamchatka, Alaska-Aleutians, Sumatra-Java, South American, and Cascadia subduction zones, which are known to produce Mw >9.0 earthquakes (38). (ii) The trench-parallel gravity anomaly (TPGA) in the North Luzon Trough (latitudes 18° to 21°) is characterized by strongly negative TPGA values of about −100 milli-Galileos (15), and previous studies suggest that great subduction earthquakes occur predominantly in areas that have a strongly negative TPGA, whereas regions with strongly positive TPGA are relatively aseismic (39). Considering all these factors together, we have no reason to rule out the possibility that zone III could generate Mw 9.0 earthquakes. There is also scant geological literature (4042) to support the inference of past basin-wide tsunamis (summarized in fig. S1).

Stochastic earthquake source model

For the fault geometry, we used the one constructed by Hsu et al. (15), which is based on earthquake focal mechanisms and the curvature of the Manila Trench (Fig. 1A). The three-dimensional fault geometry is 1200 km in length and 100 to 160 km in width, constrained by limiting the rupture depth to shallower than 50 km (50-km depth covers more than 95% of the thrust-type events capable of generating appreciable tsunami waves). The variation of fault width is mainly caused by the changing dip angle from north to south. In offshore northern Luzon, the dip angle varies from 10° to 20° at shallow depths to around 30° at depths larger than 30 km. In the southern portion of this geometry, the dip angle steepens from 10° in zone II to 30° in zone III at depths shallower than 50 km at latitude 15°N. Uncertainties do exist for the fault geometry; however, a sensitivity test of different fault geometries suggests that the tsunami impact in south China is not sensitive to the fault geometries (6). The current fault geometry and the other three constructed on the basis of previously published studies (11), earthquake focal mechanisms, and an ideal planar assumption with a constant dip angle of 17.9° produce almost identical peaks nearshore tsunami amplitudes with different return periods.

For each synthetic earthquake with a given magnitude, we followed the same approach used by Li et al. (6) to generate heterogeneous slip distribution and randomly chose a rupture area from a predefined rupture area database. A detailed description of the approach is provided in (6); here, we only summarized the key steps briefly.

The rupture dimensions (length and width) of a given magnitude were calculated by scaling relations (43, 44) and approximated to values that are multiples of 10 km. The rupture area contains a specified number of patches in the strike direction and dip direction. For example, a Mw 8.2 earthquake has 20 patches in the strike direction and 8 in the dip direction, and a Mw 8.4 earthquake has 30 and 10 in the strike and dip direction, respectively.

To generate the rupture area database, we first divided the Manila Trench area into small patches with a size of 10 km by 10 km and numbered these patches following the sequence of the strike direction first and then the dip direction. With the sequentially numbered patches, we grouped the patch numbers for each earthquake magnitude based on the required number of patches and then numbered the groups. The locations of the rupture area were randomly selected from the group numbers with the assumption that each area had equal rupture possibility.

For the heterogeneous slip distribution, we first used the hybrid model (45) (access the code at to generate high-resolution heterogeneous slip distributions in a rectangular area (the length and width of the area are equal to the rupture area calculated by scaling relations for each earthquake magnitude) and then projected the generated slip distribution to the chosen rupture area with coarser patch sizes by averaging the high-resolution slip distribution and scaling the slip value to match the total overall slip imposed by the given moment magnitude. Figure S2 (A to F) gives examples of slip distributions of Mw 8.0 to 9.0 earthquakes produced by the aforementioned method. Using this method, we generated all the 1886 earthquake events with varied slip distribution patterns and a randomly chosen rupture area.

Tsunami simulation

We used COMCOT (46) to simulate wave generation [by Okada’s model (47)], propagation, and inundation. Linear shallow-water equations were solved in deep oceans in a spherical coordinate system, and the effects of the Coriolis force were included; the nonlinear shallow-water equations were used to simulate tsunami wave propagation nearshore and inundation onshore (48). The setup of computational grids is shown in Fig. 1. Four nested grids (grids 01 to 04) were used, with the grid resolution varying from ~2000 m in the source region to 25 m in the city area. In the deep-ocean area, we used a 30–arc sec grid (ca. 925 m) from the General Bathymetric Chart of the Oceans (GEBCO) digital bathymetry/topography data set to create the computational grids for grids 01 and 02. Grid 03 has 100-m resolution covering the Pearl River Delta, and grid 04 has 25-m resolution covering the Macau area. The topographic data were derived from the following sources: (i) high-resolution bare-ground topographic data for the Macau Peninsula, the Taipa Island, the Coloane Island, and the Cotai strip measured by the Macau Cartography and Cadastre Bureau; (ii) 1–arc sec Shuttle Radar Topography Mission data for covering the Pearl River Delta; (iii) and 5-m elevation specified for the two artificial islands A and B, which are still under construction. The bathymetric data were derived from the following sources: (i) 36 nautical charts covering the Pearl River Estuary with scales ranging from 1:5000 to 1:250,000 and (ii) 30–arc sec GEBCO data (The GEBCO_08 Grid, available at, accessed in October 2014). In the innermost layer, we specified a constant Manning’s roughness of 0.025 for the inland area and 0.013 for the offshore area. The Manning’s coefficient of 0.025 was chosen after several previous studies [for example, (49)] in which coefficients ranging from 0.02 to 0.03 were used to model the tsunami inundation for coastal cities in New Zealand and the U.S. East Coast.

The tide level during a tsunami can have a major impact on tsunami inundation. Macau has a mixed semidiurnal tidal cycle in which two high and two low tides of different sizes appear every day. The maximum tidal range observed in Macau is 2.86 m, while the difference between MHW and mean low water (MLW) is 1.07 m (estimated from the tidal records during 2000–2012). Here, we followed the common practice in preparing the probabilistic tsunami inundation maps (50) and used MHW (= MSL + 0.53 m) as the tidal stage for all inundation simulation. For the two SLR conditions, we assumed that the tidal range was the same with current sea level; hence, the MHWs were also shifted 0.53 m above the adjusted sea levels.

The standards of Macau’s coastal defenses vary significantly primarily as a result of changing land reclamation practices through time. Newly reclaimed lands commonly have higher defense standards. A major portion of Macau’s coastline is protected by increasing the land elevation itself instead of using elevated vertical seawalls, which could limit the public access. The newly reclaimed lands generally have higher land elevation [>4 m above present mean sea level (PMSL)]; for example, all land reclaimed after 2001 in the airport, the east coast of Cotai, and the southern coast of Macau Peninsula (see the elevation map of Macau in fig. S8 and Fig. 1C for the time of reclamation) lies more than 4 m above PMSL. The two-dimensional simulations address the hydrodynamic characteristics of the inundation as there is no obvious “overtopping” process. On the other hand, the “older” land has lower defense standards (commonly <2 m with narrow vertical concrete walls). For example, the inner harbor area, which is the most vulnerable spot in Macau, is only protected by discontinuous low walls (with gaps in between). This area is frequently flooded by storm surges when typhoons strike Macau. The discontinuous low walls likely have limited impact on the simulation results.

Statistical analysis

For each event, we obtained the maximum inundation depths at all inland computational cells (25-m grids) by subtracting the Digital Elevation Model (DEM) data from the simulated maximum wave heights. From the thousands of tsunami inundation depth maps of Macau at different sea levels, we performed all types of statistical analysis (for example, mean and percentile of the selected group of data). The annual exceedance rate for inundation depth is calculated for each inland point with values between 0 and maximum water depth using 0.1-m increments. The water depths for certain return periods (for example, 100, 500, and 1000 years) can then be obtained through linear interpolation. Using this approach, we generated inundation maps with different return periods. For the inundation maps produced by earthquakes with certain magnitudes, we calculated the 10th, 50th, and 90th percentiles of the maximum inundation depths for each specific earthquake magnitude at each inland point. For example, for Mw 8.6 earthquakes, we first calculated the maximum inundation depths for all the 108 scenarios, and then for each inland point (one computational cell), we can get the percentiles from these 108 values. Putting all the calculated values with the same percentiles (for example, 90th percentile) together, we obtained inundation maps generated by Mw 8.6 earthquakes. The same approach was used for other sea-level conditions.


Supplementary material for this article is available at

Fig. S1. Locations of historical tsunamigenic earthquakes or tsunami-affected areas in the SCS.

Fig. S2. Examples of earthquake slip distributions.

Fig. S3. Fiftieth percentile inundation maps generated by earthquakes with different magnitudes.

Fig. S4. The variability of tsunami impact caused by Mw 8.6 and 8.8 earthquakes.

Fig. S5. The return period ratio of the six representative locations experiencing certain inundation depths at different sea levels.

Fig. S6. Tsunami inundation maps for different earthquake frequencies.

Fig. S7. Worst-case scenario.

Fig. S8. Tsunami inundation maps for 1.5-m SLR.

Fig. S9. Land elevation of Macau.

Table S1. Return periods of certain inundation depths at the six selected locations and the inundation depths with a 1000-year return period.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We appreciate K. M. Mok from the University of Macau for sharing the high-resolution topographic data of Macau and H. K. In for the tide data analysis. We thank our colleagues D. Lallemant and K. Morgan for constructive comments on this manuscript. Funding: This research was supported by an AXA Research Fund Post-Doctoral Fellowship under the project “Probabilistic assessment of multiple coastal flooding hazards in the South China Sea under changing climate” to L.L. A.D.S. was supported by the MOE Academic Research Fund (AcRF) Complexity Tier 1 Project RGC4/14 “Preparing Asian mega cities for changing climate and the potential increase in extreme sea levels and storm surges.” R.W. was partially supported by the NSF under grant no. NSF-GLD-1630099. This work is a contribution to the project IGCP 639 “Sea-level change from minutes to millennia.” Author contributions: L.L., A.D.S., and Y.W. designed the study. L.L. conducted the study and wrote the manuscript with contributions from all other co-authors. All authors contributed to the discussion and interpretation of the results. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The navigational charts in the Pearl River Estuary were purchased from Beijing Situo Ocean Information Technology Co. Ltd. Readers who are interested in the charts may contact and purchase the charts from Beijing Situo Ocean Information Technology Co. Ltd. The GEBCO data used in this study were downloaded from in October 2014. All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article