Delft3D is a process-based model that combines a short wave model (SWAN), a depth-averaged flow module, a sediment transport module (van Rijn, 1993), and a bed level update scheme based on the 2D sediment continuity equation. Delft3D has an in-built morphodynamic accelerator MORFAC (Roelvink, 2006; Ranasinghe et al., 2011) to enable faster computations. The MORFAC approach represents the lag between hydrodynamic forcing and morphological response, by multiplying the bed levels computed after each hydrodynamic time step by a time varying or constant factor (MORFAC), which enables much faster computation. The upscaled new bathymetry computed after each hydrodynamic time step is used for the subsequent hydrodynamic time step. The model structure is shown in Figure 1A. The main output from Delft3D utilized in this study is the simulated alongshore variable mean water line position along the coastline of the study areas.

The scale-aggregated model SMIC (Ranasinghe et al., 2013) provides rapid 100-year projections of coastline change adjacent to small inlet-interrupted coasts due to the combination of climate change driven sea level rise (SLR) (including Bruun effect and SLR driven basin infilling) and climate change driven variations in rainfall/runoff (including basin volume change and fluvial sediment supply effect). The schematic diagram of the model is shown in Figure 1B. The total “potential” coastline change due to the four main physical processes is calculated using Equation 1.

where, △C_T is the total potential coastline recession (m), △C_BE is the coastline change due to the SLR driven Bruun effect (m), △C_BI is the coastline change due to the SLR driven basin infilling effect (m), △C_BV is the coastline change due to the riverflow driven basin volume change effect (m), and DC_FS is the coastline change due to the rainfall/runoff driven fluvial sediment supply change effect (m). The main output from SMIC utilized in this study is the projected alongshore average coastline recession by 2100 at the study areas.

This study focusses on two real-world STIs located along the South West (SW) coast of Sri Lanka: Negombo lagoon – permanently open, locationally stable inlets (Type 1) and Kalutara lagoon – permanently open, alongshore migrating inlets (Type 2). Type 3 STIs (seasonally/intermittently open, locationally stable inlets) were excluded from this study as the available modelling approaches do not allow the quantification of coastline recession and alongshore variability of inlet-adjacent coasts for systems that are only open for some part of the year, as is the case with Type 3 STIs. This study adopts the same schematised bathymetries developed for the two study areas in Duong et al. (2017).

The two study areas on which the schematised bathymetries are based on are shown in Figure 2. In general, Sri Lanka has a monsoonal climate with a Northeast monsoon and Southwest monsoon operating respectively in the months November–February and May–September. The country’s coastal environment is classified as micro-tidal, wave dominated with a mean tidal range of about 0.5 m and an average offshore significant wave height of 1.12 m (Duong et al., 2018). The SW coast, where the two study areas are located, experiences energetic wave conditions during the SW monsoon (significant wave heights of 1–2 m, incident from the SW-W octant). Almost the entire coastline of Sri Lanka contains sandy beaches with grain diameters (D₅₀) varying between 0.2 and 0.45 mm depending on the location (Jayathilaka, 2015).

The first step in the schematised modelling approach is to generate a schematic, simplified bathymetry that follows the main characteristics of the target real-world system. The key system dimensions determined from available literature and Google Earth, and used in creating the schematised bathymetries (in Figure 3), are shown in Table 1. The schematised STI bathymetries consist of a rectangular flat-bed estuary/lagoon connected to the ocean through a straight channel. The near-shore zone consists of shore parallel depth contours, down to 10 m depth, following a Dean’s equilibrium profile corresponding with the local median sediment diameter of the area (also given in Table 1). Riverflow was introduced to the estuary/lagoon based on observation at the target STI as a flow discharge at the landward boundary of the domain. Model grid resolutions in both cross-shore and alongshore directions were optimised for Delft3D computations. High resolution (∼ 10 m × 10 m) grid cells were used in the (approximate) surf zone and inlet channel, to ensure that key physical processes in the vicinity of the inlet entrance and channel were accurately resolved by the model. The flow computational domain was extended 5 km either side of the inlet to avoid boundary instabilities affecting the area of interest (i.e., the inlet entrance), and the wave domain was created larger than flow domain to avoid any wave shadowing effect at lateral boundaries.

The schematised process-based modelling approach used here using Delft3D has been successfully applied in various previous studies such as: Marciano (2005), Dastgheib et al. (2008), van der Wegen and Roelvink (2008), Bruneau et al. (2011), Dissanayake et al. (2012), Nahon et al. (2012), van Maanen et al. (2013), Zhou et al. (2014), Nienhuis et al. (2016). The implementation of this approach to investigate climate change effects on STIs is described in detail in Duong et al. (2017) and therefore only a summary is provided here for the sake of completeness. Essentially, this approach entails the application of Delft3D to flat bed, schematised estuary/lagoon bathymetries that are constructed based on published literature and other freely available sources (e.g., aerial photographs and Google Earth), and with present and future forcing also obtained from freely available sources.

A Delft3D Present Simulation (PS) was undertaken for each schematised system, with contemporary harmonic tides together with monthly average wave and riverflow forcing (derived from measurements – see Duong et al., 2017) to reproduce the general contemporary morphodynamic behaviour of the system, to gain confidence in the model’s ability to simulate system morphodynamics. In all Delft3D simulations employed here, a MORFAC of 30 was used to capture the annual cycle of hydrodynamic forcing within a 1 year long morphodynamic simulation. The PS model has been qualitatively validated against available aerial/satellite images of the study area and compared with empirical relationships including A–P relationship (O’Brien, 1931, 1969; Jarrett, 1976), Escoffier curve (Escoffier, 1940), and Bruun inlet stability criterion (Bruun, 1978) in Duong et al. (2017) and is therefore not further described here, however, for the sake of convenience, key validation results presented by Duong et al. (2017) for the two study sites are summarised in the Supplementary Material. The validated model was implemented in a series of future snap-shot (1-year long) Delft3D simulations (CS) to investigate climate change impacts on the system (at year 2100). Here, the future forcing used represents a worse-case climate scenario following RCP 8.5, and climate change driven variations in mean water level (MWL) (i.e., SLR) and other forcing (i.e., wave conditions and riverflows) were obtained from freely available, coarse resolution sources [such as: (Stocker et al., 2013a,b), global projections by Hemer et al. (2013); see Table 2]. The basin infilling effect due to SLR was also taken into account in an offline fashion in CS simulations that included SLR. All PS and CS simulations were undertaken for a 1 year duration to represent the annual cycle of forcing (including riverflows and wave conditions).

Two sets of CS simulations were undertaken for each system to investigate: the first set investigated the impact of climate change driven variations in individual system forcings (i.e., SLR, wave height, wave direction, and riverflow), and the second set investigated the impact of climate change modified key physical processes [e.g., LST – M, ebb tidal prism (including riverflow effects) – P]. In total, 15 CS simulations were undertaken for each of the two study systems.

Although Delft3D (and other similar coastal area models) has good skill at simulating bed level changes below MWL, it is not designed specifically to simulate changes in the MWL contour (i.e., shoreline) itself, and thus the Delft3D computed coastline (i.e., the zero meter depth contour) positions are not very reliable. Therefore, the momentary coastline (i.e., MKL) approach (or in Dutch, Momentary Kustlijn – MKL) philosophy (van Koningsveld and Mulder, 2004) was adopted here to compute the alongshore variability of the coastline from Delft3D output. The MKL approach defines the coastline position as a function of the sand volume in the near shore zone, which is calculated along individual cross-shore profiles based on the volume of sand per unit alongshore length between two pre-defined, stable horizontal planes that are located above and below mean low water (MLW) and bounded by a fixed vertical landward boundary, as shown in Figure 4. In each simulation undertaken here, the standard deviation of the MKL positions calculated at all the individual cross-shore profiles was taken as an indicator of the alongshore variability in coastline position under the conditions simulated.

In this study, the momentary coastlines MKL were calculated for each STI Type at the end of each PS and CS Delft3D simulations (i.e., at the end of 1 year simulations). The (spatial) standard deviation of this end-MKL is then taken to be representative of the alongshore variability in coastline position.

As Delft3D is not an ideal tool to accurately compute alongshore average coastline recession, the reduced complexity model SMIC was used to obtain estimated of alongshore average coastline recession (by 2100, relative to the present) due to the combined effect of SLR driven Bruun effect and basin infilling, changes in riverflow volume and fluvial sediment supply (note: SMIC does not provide any information on the alongshore variability of the coastline). The model was applied to the two STIs using values obtained from published literature as shown in Table 1.

The standard deviations of the Delft3D derived MKLs were then plotted around the SMIC predicted alongshore average future coastline positions to illustrate alongshore variability in coastline position resulting from the considered climate change driven variations in system forcing.

Future mean sea level, wave conditions and riverflow at both sites are obtained from global scale projections, which by necessity are of relatively coarse scale. As such these projections may not fully represent regional features of the climate such as monsoons and regional variations in SLR (Duong et al., 2016; Ranasinghe, 2016).

The present and future alongshore variability of the coastline is here determined using Duong et al.’s (2017) process-based snap shot modelling approach. In this approach year-long simulations are performed at the two ends of the study period, as a full 100 years Delft3D simulations with concurrent and realistic tide, wave and riverflow forcing is still not possible (Ranasinghe, 2016, 2020). In this approach, while the bed level changes that may occur within the estuary/lagoon between the PS and CS are represented in an aggregated and offline way (see Duong et al., 2017), coastline plan shape changes that may occur between the PS and CS are not represented. Such plan shape changes may occur, for e.g., due to beach rotation (Ranasinghe et al., 2004; Harley et al., 2011), evolution and growth of local singularities (Ashton et al., 2001) that may result from extreme events such as very low return period storm surges (new inlet formation, extreme localized shoreline retreat) and/or riverflow/fluvial sediment supply events (formation of a new ebb delta).

Being a reduced complexity model, SMIC represents the main physical processes governing STI behaviours using aggregated equations that consider the estuary/lagoon and the inlet-adjacent coast as connected but discrete elements of an STI system. As such SMIC does not simulate local singularities such as bars, shoals, channels, shoreline protrusions/indentations that may occur and grow (e.g., due to positive feedback loops) over the 100 years simulation period. Furthermore, as in the Delft3D snap shot approach, the SMIC application here also does not account for possible long-term signatures of extreme events on mean shoreline position. SMIC also does not account for possible effects any spatio-temporal variations in LST on alongshore average coastline recession.

The alongshore average coastline recession projections obtained from SMIC are compared here with those presented by Vousdoukas et al. (2020). However, there are several noteworthy differences between the two modelling approaches. For example, SMIC takes into account the effect of estuary/lagoon processes such as basin infilling and fluvial sediment supply would have on coastline recession, which Vousdoukas et al. (2020) does not. Furthermore, Vousdoukas et al. (2020) assumes that ambient shoreline change trends would continue as they are today during the 21^st century, which SMIC does not. Thus, while the effects of processes such as alongshore gradients in longshore transport on coastline recession would be included in Vousdoukas et al.’s (2020) projection, they are not accounted for in SMIC. Finally, the SMIC application here considers a single SLR value of 1 m by 2100 [representing the high end of the likely range of global mean SLR given in Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report (AR5) (Stocker et al., 2013b) and Sixth Assessment Report (AR6) (Masson-Delmotte et al., 2021)], while Vousdoukas et al. (2020) treats SLR in a more sophisticated way using probability density functions (PDFs) of SLR for RCP 8.5 (as well as RCP 4.5) and also accounting for regional variations in SLR within a Monte Carlo simulation.