Assessing the impacts of sea level rise on salinity intrusion in the Kelani River, Sri Lanka

The global mean sea level is rising at an alarming rate of approximately 3.2 mm/year, a trend that has accelerated over the past two centuries (Church et al. 2013). The Intergovernmental Panel on Climate Change (IPCC) emphasizes that this increase poses significant risks to coastal regions, where about 10% of the world's population lives within 10 m of current sea levels (Carrasco et al. 2016). Consequences of sea level rise (SLR) include increased flooding, intensified erosion, and greater salinity intrusion into freshwater systems (Nicholls et al. 2011). These changes can severely impact estuarine water quality, potentially rendering it unsuitable for drinking, agricultural, and industrial uses (Renaud et al. 2015).

Saltwater intrusion is driven by the density gradient between saltwater and freshwater, facilitating estuarine circulation and promoting the landward movement of saline water (Hoagland et al. 2020). Several studies indicate that SLR will exacerbate saltwater intrusion in estuaries (Hilton et al. 2008; Bhuiyan & Dutta 2012), also leading to changes in stratification (Hong & Shen 2012). Therefore, determining salinity distribution along coastal rivers and predicting future distributions with respect to SLR is a major concern for water managers in estuaries (Chithra et al. 2021).

Various methodologies have been employed to analyze salt wedge transport in estuaries. While analytical models can address some limitations of field measurements, they often struggle to accurately define eddy viscosities and diffusivities (Prandle 1985). Consequently, numerical models have emerged as effective tools for analyzing saltwater intrusion (Veerapaga et al. 2019). Computer-based numerical models are well regarded for simulating complex hydrodynamic (HD) processes (Bhuiyan & Dutta 2012).

Most existing research has utilized one-dimensional advection-dispersion models due to data constraints. For instance, Bhuiyan & Dutta (2012) applied such models to the Gorai River network to assess the impact of SLR and found that a rise of 59 cm in sea level resulted in an increase of 1.5 ppt/m of SLR at a distance of 80 km upstream. Similarly, Hull & Tortoriello (1979) used a one-dimensional model to evaluate the impact of SLR in Delaware Bay, finding that salinity increased by 0.4 Practical Salinity Units (PSU) with a SLR of 0.13 m. More advanced two-dimensional models, such as Mike 21 (Mohal et al. 2006) and Delft3D (Sumaiya et al. 2015), have also been developed. However, some studies note limitations in capturing the dynamic behaviour of saline plumes due to the absence of a three-dimensional domain (Lee et al. 2016). Rice et al. (2012) developed a three-dimensional salinity model for the James River and found that the saline plume intruded 10 km further with a mean SLR of 1.0 m. Weng-Cheng and Hong-Ming also developed a three-dimensional salinity model to study the impact of SLR on salinity intrusion and transport time scales in Taiwan's estuaries.

Salinity intrusion dynamics in estuaries are closely linked to morphologic, topographic, HD, and tidal characteristics (Sumaiya et al. 2015). The interplay of tidal fluctuations and freshwater inflow is critical; upstream water withdrawals are reducing river discharge while rising sea levels facilitate further salinity intrusion (Akter et al. 2016). Numerous studies have reported the influence of SLR on salinity, residence time, material transport processes, and other relevant factors in coastal estuaries (Bhuiyan & Dutta 2012; Hong & Shen 2012). Studies conducted on the Kelani River in 1993 indicated that unacceptable salinity levels would occur if river discharge downstream of Ambatale fell below 33 m³/s, with saline wedges extending approximately 15 km upstream and fluctuating by around 4 km due to tidal effects (Nanseer & Rajkumar 2006). During the dry season (January–March), bottom salinity levels often exceed 25 PSU (ppt) (Rammadugala et al. 2007). Furthermore, a recent study by Siriwardena et al. (2021) found that the effect of backwater extends about 10 km upstream from the constructed temporary submerged barrier near the intake area.

However, no studies have examined the variation in the spatial and temporal distribution of salinity levels in a stratified domain within a river or the impact of SLR on these variations. Additionally, no salinity models have been developed for data-scarce regions, particularly for three-dimensional domains. Furthermore, no studies have been conducted to investigate how river flushing discharges vary along a river during the dry period in coastal rivers. The objective of this study is to investigate the variation in the spatial and temporal distribution of saline plumes in a data-scarce river during the dry season and to predict their behaviour under SLR scenarios using a well-developed HD and salinity transport model in a three-dimensional domain. The model will be validated using water level, discharge, and conductivity data to ensure its accuracy and reliability. The model was then applied to the Kelani Lower Basin to calculate the salinity distribution, and the results were used to assess how SLR affects salinity distribution in this domain. The impact of climate change was evaluated by considering worst-case SLR scenarios projected by the IPCC, utilizing Representative Concentration Pathways (RCP) scenarios (Katupotha 2018), and comparing these predictions with historical data trends from 1975 to 2005 (Ministry of Environment, Sri Lanka 2022). Furthermore, a practical method has been proposed to determine the tidal penetration of coastal rivers, which have highly complex dynamics, and to assess the spatial variation of flushing discharges along the river stretch during the dry period.

The structure of this paper is as follows: In Section 2, the methodology – including the study area, data collection, boundary conditions, model descriptions, model setup, simulation time period, model calibration and validation, and model results – is presented in Section 3. The discussion of the results is provided in Section 4, and the conclusions are in Section 5.

Sri Lanka is an island country in South Asia, located well within the tropics, between the northern latitude of 5° 55′ and 9° 51′ and the longitude of 79° 41′ to 81° 53′, which lies in the Indian Ocean, southwest of the Bay of Bengal, separated from the Indian Peninsula (Silva et al. 2016). Rainfall in Sri Lanka varies significantly due to the development of extreme low-pressure conditions in the Bay of Bengal. The climate on the island can be characterized into two distinct seasons: Northeast Monsoon season (December–March) and Southwest Monsoon season (May–September) (Rajapaksha et al. 2020). During the droughts of each year, when the mean sea level is high during the period from January to March of any given year, the saltwater plume intrudes more than 15 km from the river mouth near the Ambatale intake area (Ranmadugala et al. 2007). Taking into account the area most affected by the saline plume in the Kelani River due to the lowering of the riverbed especially as a result of excessive sand mining during the dry season (Siriwardena et al. 2021), the Lower Kelani River Basin has been chosen for this study, bounded by Hanwella and the Modara sea outfall. Lower boundary conditions were set at the sea outfall, incorporating hourly tidal fluctuations (represented as a stage) (Samarasinghe et al. 2022). Consequently, the stage hydrograph at Nagalagama Street successfully replicated tidal fluctuations, including the backwater effect (Bakhtyar et al. 2020). The upper boundary was established by considering locations where no salinity plume was observed and where existing river gauge stations are located, particularly in the Hanwella area.

The data required for setting up the model was obtained through collaboration with relevant authorities and departments such as the National Water Supply & Drainage Board, the Irrigation Department, the Survey Department, and the Meteorological Department. Time series data were collected to enhance the model's precision. This dataset included a variety of information such as topography (digital elevation model – DEM), canal cross-section data, river water levels and discharges, details about existing water bodies, tidal water levels, temperature records, salinity measurements, humidity levels, wind speed, and wind direction. Descriptive statistics have to be analysed of the measured data samples, as shown in Table 1, which clearly indicates central tendency and dispersion, along with an assessment of the data quality related to the data used in this study (Burgan et al. 2013).

A detailed description of the model setup, calibration, and approach used can be found in the research paper by Melo et al. (2020). The numerical model selected was the Delft3D software, which is widely used for simulating HD processes in estuaries (Iglesias et al. 2019). This model operates by solving the continuity and momentum Equations (1)–(4) for an incompressible fluid, as outlined in the research by Broomans (2003).

Delft3D is a unique and fully integrated computer software suite aimed at facilitating a multi-disciplinary approach and conducting 3D computations for coastal and river environments to execute HD models and water quality assessments (Deltares Delft3D-Flow 2013). Initially, HD conditions such as velocities, water elevations, density, vertical eddy viscosity, and vertical eddy diffusivity were computed using Delft3D-FLOW. Subsequently, the output from Delft3D-FLOW served as input for the D-Water Quality model (Deltares Delft3D-Flow 2013). Delft3D-Flow possesses the capability to calculate non-steady flow and transports resulting from tidal and meteorological forcing in both 2D and 3D domains (Deltares Delft3D-Flow 2013). DELFT3D–water quality (WAQ) can simulate a wide range of substances, solving mass balance equations encompassing advection, dispersion, and reactions (Deltares D-Water Quality 2013).

To facilitate salinity simulations, a 3D model was employed, incorporating four vertical sigma layers within the same domain utilized in the HD model. Experimentation was conducted with finer vertical resolution, yet it was determined that increasing the number of layers further had no significant impact on the numerical outcomes due to the shallowness of the domain (Matsoukis et al. 2022).

Salinity holds a unique position in D-Water Quality, as it is not subject to water quality processes but solely to transport mechanisms. Salinity was created as the substance to be run by advection-dispersion formulae in the model domain. Consequently, in the transport model domain, four equal layers were vertically aggregated using the layer editor option in the Delft3D–WAQ model, aligning with the same domain utilized in the HD model for the coupling process. The same methodology used in this study, including the coupling of water quality and HD models in 1D, 2D, or 3D, can be applied to analyse similar problems of a comparable nature.

It is evident from the data that the minimum discharge spells occur predominantly from January to April/May each year. Consequently, January and February (January 22, 2020–February 7, 2020) were selected as the simulation period for this study to encompass the worst-case scenario for the salinity model.

Model calibration and validation are essential steps in mathematical modeling, typically involving dynamic parameters such as water levels, current speeds, and salinity levels (Shukla et al. 2015). Before proceeding to verify the integrated model, it is essential to first verify the HD model independently. The performance of the HD model was evaluated by comparing its results with observed water levels at the Nagalagama Street River gauging station. Calibration parameters, including the Manning coefficient, horizontal eddy viscosity, and vertical eddy viscosity, were iteratively adjusted until the predicted water levels at Nagalagama Street closely matched with the observed values for the period from January 22, 2020 to February 7, 2020.

Once the calibration parameters for the HD model were optimized, the validation process involved comparing water level predictions for two additional periods: from January 10 to January 17, 2020 at Nagalagam Street, and from January 22 to February 7, 2020, at the Ambatale intake area and Hanwella. Subsequently, a well-tuned HD model was utilized to couple with the hybrid model, and the integrated model underwent calibration and validation using measured salinity values collected from the National Water Supply and Drainage Board.

The simulation period selected for the calibration process was from January 22, 2020 to February 7, 2020, as measured salinity levels were available during this timeframe. The comparison was conducted using the measured salinity levels in the top and bottom layers at the Kelanisiri bridge area, obtained from the National Water Supply and Drainage Board's Ambatale Water Treatment Plant (WTP) Laboratory. Model parameters for horizontal eddy diffusivity and vertical eddy diffusivity were adjusted until the predicted salinity levels were closely aligned with the measured values. This simulation and comparison were performed considering two layers within the river to assess the model's performance in a three-dimensional domain. Throughout the calibration procedure, two layers were considered, representing the river's bottom and top layers, as salinity measurements are typically taken at these two levels by National Water Supply & Drainage Board (NWSDB).

The adverse impact of Global Climate Change on saltwater intrusion processes in coastal aquifers is exacerbated by SLR (Chang et al. 2011). The extent of this rise primarily hinges on the rate of future carbon dioxide emissions. The IPCC Fifth Assessment Report (AR5) delineated projections for 21st-century SLR based on varying levels of greenhouse gas emissions. These projections were classified into four scenarios under RCPs, reflecting diverse emission trajectories: RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5. For this study, the three most severe scenarios were selected: RCP 2.6, RCP 4.5, and RCP 8.5. Additionally, two time spans were chosen to represent the 21^st century: 2050, signifying the mid-century, and 2090, symbolizing the century's end. A summary of the selected scenarios is presented in Table 2 (IPCC 2014; World bank n.d.).

To comprehensively assess the spatial variations in salinity concentration along the river, the study employed modelling techniques to estimate the average salinity concentrations across different layers spanning from the base year to projected future scenarios. The primary objective was to delineate the distances over which salinity intrusion occurred along the river. Salinity intrusion length specifically denotes the distance from the mouth of the estuary to the furthest upstream point where the mean salinity level at a given cross-section is observed.

The propagation of the resultant forces from tides and freshwater through coastal aquifers, in both landward and seaward directions, is referred to as tidal penetration. In estuaries and tidal rivers, tidal penetration refers to the extent to which tides can transport saline plumes upstream. This phenomenon is influenced by two primary factors: tidal cycles and upstream river discharges.

Studying tidal penetration lengths in HD models of rivers is challenging due to the dynamic nature of boundary conditions. To address this, a methodology was employed to identify tidal penetration in the lower Kelani basin. Initially, a constant river discharge was set at the upstream end of the model, while a temporal tidal cycle was applied at the downstream end. The model was iterated until the river pattern stabilized, indicating the most landward extent of tidal penetration. This process was repeated with various constant upstream river discharges. The analysis of tidal penetration in a river can thus determine how tidal penetration lengths vary with different upstream river discharges.

Furthermore, in addition to visual observation, sensitivity analysis of the results can be conducted using statistical indicators. The indicators employed include root mean square error (RMSE), goodness-of-fit (R²), percentage of deviation from observed results (PBIAS), and Nash–Sutcliffe efficiency (NSE) (Bhuiyan & Dutta 2012; Edirisooriya et al. 2022).

Optimal model performance is anticipated when RMSE values approach zero and R² values approach unity (Bhuiyan & Dutta 2012). PBIAS assesses the average tendency, where positive values signify model underestimation bias and negative values indicate model overestimation (Gupta et al. 1999). NSE (Nash & Sutcliffe 1970) is a normalized statistic that gauges the relative magnitude of residual variance in comparison to measured data variance. NSE ranges from −∞ to 1 (inclusive), with an NSE value of 1 representing the highest efficiency of the model. Efficiency closer to 1 indicates greater model accuracy, while values between 0 and 1 generally denote acceptable performance levels. Values ≤ 0 indicate unacceptable performance (Nash & Sutcliffe 1970). Considering more than the visual observation, the statistical evaluation also confirms the good establishment of the integrated model, as mentioned in Table 3.

A well-developed HD c and advection-dispersion three-dimensional model was employed to simulate the behaviour of the saline plume in the lower Kelani basin across six established scenarios, relative to the present condition in 2020, particularly during dry spells each year. It was observed that a higher concentration of salinity intruded at a greater percentage during spring tides compared to neap tides. The intrusion scenario primarily depends on the aforementioned tidal penetration and Savenije's analytical equation (Savenije 2001), which was derived from the complete St. Venant equations in a Lagrangian reference frame. This equation accurately describes how river discharge has a negligible influence on tidal damping in the downstream sections of river estuaries compared to its influence upstream, as the cross-sectional area is larger downstream than upstream, where the freshwater velocity exceeds the tidal velocity (Horrevoets et al. 2004).

It is apparent that with deteriorating environmental conditions, there is a discernible rise in predicted salinity intrusion. These predictions have been consolidated and quantified in Table 4, enhancing the depth of the analysis.

The worst-case scenario is represented by RCP 8.5, while RCP 4.5 depicts an intermediate emission pathway. Under these conditions, the average bottom salinity intrusion is projected to increase by 4.47% in 2050 and by 17.35% in 2090. For the surface layer, salinity is expected to increase by 34.58% in 2050 and by 45.45% in 2090, indicating the most severe changes. Additionally, significant rises in salinity values are projected under RCP 8.5, with increases ranging from 13.17 to 27.49% for the bottom layer and from 6.06 to 75.76% for the surface layer. These changes are expected to occur between 2050 and 2090 in the Kelanisiri bridge area.

Salinity intrusion lengths are categorized into three main types based on different tidal conditions: intrusion length during low water slack, intrusion length during high-water slack, and tidal average intrusion length, which represents an average encompassing both low and high-water slack periods (Savenije 1993; Prandle 2004). These metrics provide crucial insights into how salinity levels vary spatially along the river under varying hydrological conditions, informing strategies for managing and mitigating salinity impacts on aquatic ecosystems and water resources.

The results show a significant increase in salinity intrusion, particularly under the RCP 8.5 scenario for the year 2050 and all RCP scenarios for the year 2090. When comparing current conditions to projected future conditions, the intrusion lengths for RCP 8.5 in 2050 increased by 500 m. For RCP 2.6 in 2090, the intrusion length increased by 3.5 km. Under the RCP 4.5 scenario in 2090, the intrusion length increased by 8.5 km, while for RCP 8.5 in 2090, the increase was 11 km. These intrusion lengths are summarized in Table 6. The comparisons were made with respect to a 1 ppt threshold value of the plume.

These significant increases in intrusion length demonstrate the considerable impact of SLR on salinity intrusion in the Kelani River in the future. The findings underscore the importance of accounting for climate change scenarios in managing and mitigating the effects of salinity intrusion on freshwater resources.

The developed hybrid model analyzes the impact of SLR on salinity intrusion, specifically focusing on spatial and temporal variation in salinity plumes in the bottom and surface layers throughout the lower basin. The model findings reveal that, compared to current conditions, the average increase in bottom layer salinity is projected to be 11.8, 4.5, and 13.1% in 2050 and 13.8, 17.4, and 27.5% in 2090 under the environmental scenarios RCP 2.6, RCP 4.5, and RCP 8.5, respectively. Additionally, the simulated results indicate that a saline front of 1 ppt, which significantly impacts water treatment processes, is expected to intrude further inland by approximately 0.5, 3.5, 8.5, and 11 km for respective SLR values of 21, 32, 37, and 48 cm, relative to conditions observed in 2020.

The research clearly addresses how the salinity plume has varied along both the bottom and top layers and how the length of saline intrusion has changed along the river in relation to future sea level rise. It also examines how flushing discharge has varied along the river stretch. A significant finding of this study is that the minimum required freshwater discharge to flush out salinity near the Ambatale intake area has increased to 90 m³/s in the Hanwella area, compared to the previously determined value of 33 m³/s in 1992. This increase can be attributed to variations in tidal patterns, sea level rise, changes in riverbed levels, alterations in freshwater discharge patterns, adjustments in river cross-sections, and variations in river bottom friction values. This model has been developed in a data-scarce region, and by calibrating and validating the HD model with water levels, the coupled 3D hybrid model has been calibrated and validated using measured conductivity values from the bottom and top layers, without considering peak area predictions. The comparisons of 3D salinity were further supported by statistical analyses in addition to visual observations.

This model can be utilized by experts for HD analyses and other water quality assessments in the study area. However, the situation may worsen due to future reductions in riverbed levels and changes in upstream river discharge patterns, for which this study has limited data.

At present, the model is constrained to four vertical layers due to data availability, but increasing this number would improve accuracy. While this study focuses on the dry season, when salinity is most critical, future research could extend to include the entire year, incorporating hydrological effects from the wet season. The conductivity data used for boundary conditions assumes stability in the future; however, predicting changes in these conditions could yield more precise results.

In this study, a three-dimensional HD model has been independently developed and coupled with a salinity transport model. This HD model can be integrated with any solute transport model within the domain to predict solute behaviour effectively. Unlike previous studies, this research examines HD aspects in the model domain while also analyzing the variance of the saline wedge, providing a more comprehensive understanding of the system.

This study integrated the Delft3D-Flow HD model and the salt transport model of Delft3D–WAQ for the Kelani lower basin to simulate the behaviour of the saline plume in a three-dimensional domain, both the top and bottom layer-wise, and to analyse the variation of these factors with respect to future sea level rise. The integrated model was thoroughly calibrated and verified using observational data on water levels and salinity distributions from 2020. The simulation results showed good agreement with observed data during both calibration and validation periods, as confirmed by visual comparison and statistical evaluation.

Importantly, depth-averaged or one-dimensional water quality and hybrid models are inadequate for effectively studying saline plumes in rivers. It is essential to develop a three-dimensional domain to accurately model and understand the behaviour of saline plumes. Therefore, this study utilized the integrated model to analyse the plume's behaviour in 2020, as well as its projected behaviour for the mid-21st-century (2050) and the end of the 21st-century (2090) under three different environmental scenarios.

In this finding, it can be observed that the saline front with a density of 1 ppt is expected to intrude by 8.5 and 11 km, respectively, compared to the current extent of the 1 ppt plume, with predicted sea level rises of 37 and 48 cm. This clearly indicates that the saline plume has intruded towards the land side, which could have further crucial implications, such as riverbed lowering and changes in the upstream river patterns.

The findings clearly indicate that projected future conditions will exceed the minimum threshold salinity levels in raw water intake for the Ambatale area. Moreover, existing mitigations appear insufficient to ensure the provision of safe raw water to the treatment plant at Ambatale. Therefore, identifying and implementing a permanent solution to address this issue in the near future is imperative.

Finally, models of this type can be employed to design salinity barriers for coastal rivers by combining hydrological behaviour during the wet season to mitigate the effects of upstream flooding.

The authors express their gratitude to the National Water Supply & Drainage Board, Sri Lanka, for generously providing the data necessary for conducting research simulations. Additionally, sincere appreciation is extended to the Department of Meteorology, Survey Department, Mahaweli Authority, NARA Institute, and the Irrigation Department, Sri Lanka, for providing their recorded data throughout the study period.

V.D.A. conceptualized and validated the work, developed the methodology and software, rendered support in formal analysis, and data curation, wrote the original draft, and investigated and visualized the project. P.N. conceptualized the work, developed the resources, wrote and reviewed and edited the article, and supervised the process. K.P.P.P. conceptualized the work, developed the resources, wrote and reviewed and edited the work, and supervised the process.

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.