Monte-Carlo simulations of the impact of climate change on sea-level change are an elegant tool to incorporate projection uncertainties in a physically meaningful way. Such models produce a large set of individual realizations of sea-level time series that then can be characterized using descriptive statistics at a given time. For example, Table 1 provides different percentiles for the years 2050, 2100, and 2200 for RCP 2.6 and RCP 8.5 for different Antarctic Ice-Sheet dynamics treatments.

For our study, we created a method to generate the desired number of sea-level time series, over the period from 2000 to 2100 at a 1-year increment, from the percentiles given in the Table 1 from Kopp et al. (2017). To create these sea-level time series, which we refer to as pseudo-realizations, we fit a theoretical cumulative distribution to the provided percentiles in the provided years. Calculating the derivative of these theoretical cumulative distributions provides us with the probability density function. We sample the probability density function in each provided year using the Kernel Density Estimates method. We create an array with the size of the number of pseudo-realization for the year 2000 in which the sea level is defined as zero. Taking the sampled distributions in the years 2050, 2100, and 2200, we append values to the pseudo-realizations array for the provided years by randomly choosing a value from the sampled distribution in 2050, 2100, and 2200, respectively, for which η_sl(t_{i +1}) = η_sl (t_i ) + Δ_{i +1}η_sl where t_i = 2050,2100,2200, and Δ_iη_sl is small range of sea level. As the last step, for each of the time series containing now values for 2000, 2050, 2100, and 2200, we employ polynomial function to interpolate for the years in the between t_i .

If the number of pseudo realizations is sufficiently large (>1000), the recalculated percentiles from the pseudo realizations fall within an error of less than 2% of the percentiles provided in Kopp et al. (2017).

From a modeling viewpoint, the influence of sea-level rise and other climate-change impacts on the nearshore bathymetry and, subsequently, the tsunami hazard can be incorporated in two different ways. The first way is to raise or lower the water level or mean sea level according to RCP sea-level projections, which is an approach that is also referred to here as the bathtub approach. As mentioned above, we consider this simple adjustment as static because features, such as barrier islands or marshes, do not change their horizontal locations. In the bathtub approach, adjustments due to sea-level rise only occur vertically, where the vertical adjustment is simplified to a one-to-one shift with sea-level rise.

The second approach to incorporating sea-level rise and climate-change impacts that alter bathymetries is to employ system-level simulation tools, such as the previously mentioned LTM model by Lorenzo-Trueba and Mariotti (2017). The LTM model solves a set of ordinary differential equations (ODEs) that calculate the horizontal and vertical positions of the subcomponents boundaries of the BML system. Figure 1 depicts the subcomponents of the LTM model. The red crosses show the locations of the governing equations, while the green arrows indicate important fluxes and parameters, but also where non-linear connections between the different equations exist. The red arrows mark the location of the important parameters, and important parameters, such as the barrier height (H_b ) and width (W_b ), or mainland-marsh width (W_bm ).

Without loss of generality, the LTM model is a morphokinetic flux model, and, therefore, it does not focus on the physical, chemical, and biological processes directly to drive the development of the BML system. We refer to Lorenzo Trueba and Mariotti (2017, and references therein) for an extensive discussion about the governing equations and their inputs, initial conditions, and applications of the LTM model. As mentioned in the previous paragraph, the LTM model solves a set of coupled ODEs that describe the horizontal and vertical changes of the subcomponent boundaries (see Figure 1). The key parameters are listed in Table 2.

Yet, this set of ODEs can be solved efficiently enabling us to create a Monte-Carlo-type computational framework.

The width of the barrier island is defined as W_b = x_b − x_s , the widths of the backbarrier (W_bm ) and mainland (W_mm ) marshes with W_bm = x_bm − x_b , and W_mm = x_l − x_mm . The lagoon width is calculated by W_la = x_l −x_b .

As mentioned earlier, climate change is the main driver of sea-level rise in the LTM model. However, other climate change impacts also affect the coastal and marine physical environment, for example, the nearshore wind speed or the sediment overwash flux (for more information about both parameters, see below) and can drive significant changes in the coastal system as a whole as well as its subsystems. To gain initial insights on which parameters are the most influential, we performed a small and limited sensitivity study, varying parameters found in Lorenzo Trueba and Mariotti (2017, Tables 1, 3). The result of this preliminary testing revealed that the shoreface sediment flux, Q_fs , the overwash sediment flux, Q_ow , and the influx of offshore sediment flux, I_os are the most sensitive processes and those that have the most significant impact on, the overall dynamics of the BML system. It is interesting to note that all three parameters (Q_fs , Q_ow , I_os ) relate to climate change impacts in different ways. For example, wind speed predominantly controls the shoreface sediment flux. Eichelberger et al. (2008) and Zeng et al. (2019) study how climate change generally impacts wind speed. Tides and the concentration of nutrient-rich open-ocean sediment concentration govern the influx of nutrients into the lagoon and by fractionation the bounding marshes and impact I_os . Pickering et al. (2017) indicates that the impact of climate change on the tidal range and tidal wave period –both are a function of the shelf water depth– is at a maximum within 10% of today’s value. Therefore, the impact of future change in open-ocean sediment concentration is likely to have a much more significant impact on the influx of offshore sediments. Sediment overwash onto the barrier island and into marsh-lagoon subsystems is a necessary process. Leaving the ratio of overwash sediment flux that moves onto the barrier and into the marsh-lagoon subsystem constant, the magnitude of sediment overwash is the most significant parameter influencing the entire BML-system evolution, aside from sea-level change.

The result of our aforementioned preliminary parameter study helps us to identify the parameters that need to be studied and put into the context of future climate change. Sea-level rise is, by far, the best-studied climate-change impact. However, open-ocean sediment concentration, overwash sediment flux, and wind speed as a function of climate change are not well understood, and future predictions are not available or part of ongoing controversial discussions within the respective scientific communities. To accommodate changing values for sediment overwash flux, wind speed, and open-ocean sediment concentration, we changed the original LTM model as presented by Lorenzo-Trueba and Mariotti (2017) slightly, without changing the governing equations. Furthermore, we also adopted the value ranges of sediment overwash sediment flux and the available offshore sediment from Lorenzo-Trueba and Mariotti (2017) and wind speeds from Zeng et al. (2019).

To account for future changes and trends of the open-ocean sediment concentration, wind speed, and overwash sediment flux, we assumed that the mean values of all three parameters over the period from 2000 to 2100 could either remain constant, increase or decrease. Again, ranges for end values for these different trends are taken from Lorenzo-Trueba and Mariotti (2017) and Zeng et al. (2019). Furthermore, we sample a Gaussian distribution of wind speed, overwash sediment flux, and available offshore sediment whose mean values are provided following the trend schemes mentioned above and standard deviations defined to be 10% of the respective mean value to create an, albeit, artificial spread of values in each year. We argue that this proposed method to mimic annual fluctuations simulates their stochastic nature as a Wiener process reasonably well.

Lastly, the geometry, especially of the barrier island, shown in Figure 1, is simplified, and the barrier does not feature a characteristic profile containing a dune. The barrier height as considered in the LTM model can be interpreted as an average height. Nevertheless, one could argue a profile containing a dune might have a significant influence on how a tsunami wave would overtop the barrier. Therefore, we assume that the dune is eroded by the tsunami wavefront as it travels over the island. Hence, we argue that the average barrier height is appropriate to describe the height of a barrier.

We employ GeoClaw, part of Clawpack, for tsunami simulations (no sediment transport). GeoClaw has been under development since 1994 (i.e., George, 2008; Berger et al., 2011; LeVeque et al., 2011). Clawpack solves hyperbolic systems of partial differential equations in one, two, and three dimensions. For flow problems in large spatial domains, depth-averaged Shallow-Water equations (SWE) are usually sufficient. However, some problems appear in solving the SWE, especially for extreme wave propagation and flow over variable bathymetries and topographies. For this reason, a special subset of Clawpack, GeoClaw, includes adaptive mesh refinement (AMR) (Berger et al., 2011). GeoClaw has been employed for a variety of problems in tsunami science (e.g., MacInnes et al., 2013; Melgar and Bock, 2015; Melgar et al., 2016). GeoClaw is verified and validated with all problems given in NOAA’s standards and procedures for tsunami inundation codes (George, 2008), with a recent application and benchmarking for storm surges (Mandli and Dawson, 2014).

We consider the earthquakes in this study to range from magnitudes M_w8.0 to M_w9.0 with 11 magnitude steps. We base the surface deformation field computations on random depth underneath the ocean floor between 5 km and 100 km while leaving the rake, dip, and strike angles constant. A total of 50 earthquakes in each magnitude step is used to ensure a stable tsunami-runup distribution. The area of the earthquake is computed using a standard empirical formula by Strasser et al. (2010). The earthquakes distance to the barrier island is kept constant at about 1970 km to ensure that the coastal environment is located in the far-field for any given earthquake.

As mentioned in the previous paragraph, the earthquake location is about 1970 km away from the BML system. The water depth is constant at 1000 m up to 1800 km from the earthquake area and then increases with a simple sloping beach where the shoreface in front of the barrier connects to the BML system. Figure 1 depicts the general BML system, while Figure 3 in Results below shows the BML system at scale (solid red line and a and b).

The barrier initially is 2 m high and 300 m wide, the back-barrier marsh is 4 km wide, the lagoon is 30 km long, and the mainland marsh initially is 2.5 km wide. These values do not represent any particular BML system by any means and are adopted from Lorenzo-Trueba and Mariotti (2017).

As discussed earlier, wind speed, overwash sediment flux, and open-ocean sediment concentration have profound impacts on the BML model, and, unlike sea-level rise, it is unclear how future changes under different RCPs influence the trends for these parameters. The table in Figure 6 contains the different assumed climate-change scenarios for these additional parameters and notes if the parameters remain constant, increase or decrease as described above. We acknowledge that many more combinations of constant, increasing and decreasing trends are possible among the three parameters, but argue that any additional combination can be realized by combining scenarios A-G with each other.

The parameter space of our study consists of the number of different RCP scenarios, two different treatments of Antarctic Ice Sheet uncertainty, the number of climate-change scenarios for additional parameters (wind speed, overwash flux, and sediment concentration), the number of realizations of the coastal-evolution model, the number of different earthquake magnitudes, and the number of different random depths for each earthquake magnitude. This results in a total of 2 RCPs (RCP 2.6 and 8.5) × 2 AIS models (K14 and K17) × 7 different climate-change scenarios (A-G) × 1000 (realizations per scenario) × 11 earthquake magnitudes × 100 earthquake depth per magnitude step, which totals 30,800,000 individual GeoClaw simulations on a modern Intel i7 (5th generation) processor, each GeoClaw simulation takes 4 minutes, which would result in a total computation time of 85 years. We, thus, make use of high-performance computing resources to execute these simulations.

The range between the 5^thand 95^th percentile of the sea-level distribution as a function years is depicted in Figure 2, while Table 3 contains the median (50^th percentile), the percentile range 1-99, 5-95, and 17-83. From Figure 2A, we see that for RCP 2.6 the range between the 5^th (dotted lines) and the 95^th (dashed lines) percentiles is increasing to similar values through the years. We note that K14 between 2100 and 2120 is slightly larger than K17. In 2100, the 5^th- 95^th percentile range is between 0.45 m to 1.00 m for RCP 2.6. There are more significant differences between K14 and K17 for RCP 8.5 (Figure 2). Around the year 2030, the 5^th-95^th percentile ranges diverge, resulting in a range between 0.43 m and 1.5m for K14 and between 1.06 m and 2.6 m for K17 in 2100.

From Table 3, we can see that the median values for RCP 2.6 K14 and RCP 8.5 K14 are similar for the years 2025 (RCP 2.6: 0.12 m; RCP 8.5: 0.13 m) and 2050 (RCP 2.6: 0.24 m; RCP 8.5: 0.30 m), but are significantly different for 2100 with 0.53 m for RCP 2.6 and 0.87 m for RCP 8.5. The median values for RCP 2.6 K17 and RCP 8.5 K17 for the year 2025 are the same, 0.11 m. In the year 2050 the values are 0.23 m for RCP 2.7 and 0.38 m for RCP 8.5. However a stark difference is apparent for the 2100 where the value for RCP 2.6 is 0.48 m and the value for RCP 8.5 is 1.61 m. Looking at the extreme values (95^th, 99^th, and 99.9^th percentiles), we see similar trends between RCP 2.6 and 8.5 for K14 and K17 as described for the median values.

Figure 3 depicts a comparison between the static (3a+c) and dynamic (3b+d) bathymetry models under sea-level rise. Both models are based on the same sea-level change time series and start with the same bathymetry (solid blue line). The simulations end at t_e = 100 years, which is indicated by the dashed blue line in both models.

Comparing the end states of the static and dynamic model with each other, we can see that the barrier height in the static case significantly decreases, while the barrier height in the dynamic case changed very little from the year 2000 to 2100 (Figures 3B, D). Due to the simple dependence of the static model to sea-level change, the water layer over the marshes and the lagoon increases in thickness at the same rate, and the barrier decreases its height. The depth of the marshes in the dynamic model changed some, but it seems as if the marshes were able to maintain more or less the same depth as the sea level was rising, but marsh edge erosion reduced the marsh width. Furthermore, the lagoon depth at t_e is 6m in the dynamic model, significantly deeper than in the static model (lagoon depth at t_e = 3.7m). This result in itself is interesting, but the more interesting part from a theoretical point of view is the seemingly very non-linear response to the sea-level driving change in the early years of the simulation. The yellow lines mark the half time of the simulation ( 1 2 t e = 50 y e a r s ) , and in the static model, the position of this yellow line is closer to the start state at time t_o (the year 2000) because the sea-level change solely governs the temporal development. It is projected that the sea level changes faster in the latter half of this century (> 50years) compared to earlier in this century. However, in the dynamic model, the yellow line is closer to the end at t_e , which indicates a very non-linear response of the lagoon depth to the relatively smaller sea-level change in the early years of the simulation. The lagoon depth for 1 2 t e < t < t e does not change much.

We employ the tsunami run-up as an indicator of the tsunami hazard, which is defined in our context as the vertical distance between the elevation of the maximum tsunami inundation and sea level. Tsunami run-up probability densities for all of the considered scenarios are provided in Figure 4 (RCP 2.6) and 5 (RCP 8.5). For both figures, the left column (a, c, and e) contains the static bathymetry results, and the right column (b, d, and f) the results for the dynamic bathymetry. Furthermore, the solid lines represented the run-up probability density for K14 and the dashed lines for K17; the different colors represent different years (orange: 2025, green: 2050, and blue: 2100).

Focusing first on RCP 2.6 (Figure 4), we can see that the general geometry of the run-up probability density in the case of a static bathymetry (a, c, and e) for earthquake magnitudes 8.6, 8.8, and 9.0 are very similar, but note that the run-up values along the x-axes change significantly from values between 0.075 m and 0.28 m for M_w8.6, from 0.91 m - 1.04 m for M_w 8.8, and 5.52 m - 6.8 m for M_w 9.0. The changes from 2025 to 2100 are consistent with sea-level changes (see Table 3), which is the only variable for static bathymetry. Furthermore, note how the modes of the run-up densities move toward larger run-up values, and are more broadly distributed across run-up values, with increasing time. For the dynamic case (Figures 4B, D, F), the value ranges for the run-up probability densities show a similar increase for the three different earthquakes, which is close to an order of magnitude increase in run-up range from M_w 8.6 to M_w 8.8 to M_w 9.0. For M_w 8.8 and 9.0 scenarios, the probability densities for the dynamic case are skewed toward larger values, and the range of the K14 and K17 curves in a given year have a much wider spread compared to the static bathymetry. For example, the run-up density for K17 in the static bathymetry of M_w 8.8 in 2100 ranges from 0.91 m to 1.04 m, while the run-up density interval for the dynamic bathymetry for the same case reaches from 0.2 m to 3.8 m, exceeding the range of sea-level change in the same year significantly. Furthermore, the modes of the run-up probability density in the dynamic case for all years and AIS treatments (K14 and K17) are consistent for the three earthquake magnitudes; with the exception of a slight spreading toward higher run-up values with time, the density shape does not change in the different years.

For RCP 8.5 (Figure 5), observations about the run-up probability densities for the static and dynamic bathymetry and different earthquakes are similar to that for RCP 2.6. The ranges of run-up values for the static case are consistent with the sea-level rise ranges in the different years (see Table 3).

While the run-up probability densities for the static bathymetry (Figures 5A, C, F, and Figures 4A, C, F) are only a function of sea-level rise, the dynamic bathymetry evolution in time additionally depends on offshore sediment availability, overwash sediment flux, and wind speed. Figure 6 depicts the run-up values for an M_w 8.9 earthquake for RCP 2.6 and 8.5 for K14 and K17 as box plots. As a reminder, the line within the box marks the median run-up value, while the height of the box is referred to as the interquartile range (25^th- 75^th percentiles). We employ the 5^th percentile to define the location of the lower whisker and the 95^th percentile to defined the location of the upper whisker. Note that we exclude the date that is either below the 5^th and above the 95^th percentile as outlier values. The letters refer to the climate change scenarios as provided in the table included in Figure 6. Focusing on RCP 2.6 first (Figures 6A, B), we can see that the median run-up values of the different climate-change impact scenarios vary between 2.00 m and 3.75 for 2025, between 2.00 m and 4.64 m for 2050, and between 2.30 m and 4.10 m for 2100. Note that scenario E defines the maximum of provided intervals and scenario D the minimum for RCP 2.6 for K14 and K17. Interestingly, the median run-up values for scenario E are larger for K14 than for K17 in 2025 and 2050, while K17 is larger in the year 2100. For the interquartile range, the behavior in 2025 and 2050 is very similar. The interquartile range in 2025 has its minimum range at 0.17 m (scenario D) and its maximum at 1.76 m (scenario E, K14). For the year 2050, the maximum interquartile range is 2.30 m (scenario E, K17) and the minimum is 0.25 m (scenario D, K14). The previously described pattern of minimum and maximum interquartile ranges does not continue in the year 2100. In 2100, the maximum interquartile range for K14 is at 2.10 m with scenario E; the maximum for K17 occurs at scenario F with 1.86 m. The lower and upper whiskers exhibit a similarly consistent pattern for both K14 and K17. The largest upper whiskers occur in scenario E and increase from 2025 to 2100 is similar ways. However, the upper whisker for scenario E is with 6.36 m for K17 larger than the 6.09 m for K14, which is consistent with the slightly larger sea-level range as depicted in Figure 2. For the lower whiskers, it is interesting to observe that the minimum values occurs at scenario F and is close to zero for K14 for the years 2025, 2050, and 2100; while for K17 the minimum value is only close to zero in the years 2025 and 2050.

Figures 6C, D depict the run-up values as box plots for the different climate change scenarios and years for RCP 8.5. While for RCP 2.6, the maximum of the medians can be found at scenario E, the maximum median values occur at scenario C for RCP 8.5 for K14 and K17 (2025: 3.91 m [K14], 2050: 4.72 m [K17], and 2100: 4.85 m [K17]). The interquartile-range and whisker pattern are consistent between K14 and K17 from 2025 to 2100.

Figure 7 depicts the run-up as function of the earthquake magnitude for RCP 2.6 (a & b) and RCP 8.5 (c & d). The run-up values for RCP 2.6 for earthquake magnitude smaller than 8.8 (Figure 7A) are small with neither the median nor the maximum exceeding one meter in any given year. For larger earthquakes, the run-up values change dramatically for larger earthquakes (M_w 8.8 - M_w 9.0, Figure 7B) where the medians increase from about 1 m for M_w 8.8 to values between 3.5 m and 4.1 m for M_w 8.9 to values between 7.0 m and 8.3 m for M_w 9.0. Notably, while the median is highest for 2100, the maximum values of the tsunami run-up for K14 and K17 for a M_w 9.0 earthquake occurs in the year 2050 (Figure 7B). Similarly, the median and maximum values of tsunami run-up remain below one meter for M_w 8.5 to M_w 8.7 for RCP 8.5 (Figure 7C), except for K17 in 2100 whose maximum is 1.1 m. Furthermore, the median changes for larger magnitudes (Figure 7D), where median values have similar values compared to RCP 2.6. It should be noted here that the median values for K17 in the year 2100 for a M_w 9.0 earthquake is slightly larger. The biggest difference in run-up values for M_w between K14 and K17 for RCP 8.5 occurs in 2100. The maximum for RCP 8.5 K14 is 15 m, while the maximum for RCP 8.5 K17 with 30 m twice as large.