Projections of salt intrusion in a mega-delta under climatic and anthropogenic stressors

Rising temperatures, rapid urbanization and soaring demand for natural resources threaten deltas worldwide and make them vulnerable to rising seas, subsidence, droughts, floods, and salt intrusion. However, climate change projections in deltas often address climate-driven stressors in isolation and disregard parallel anthropogenic processes, leading to insufficient socio-political drive. Here, using a combination of process-based numerical models that integrate both climatic and anthropogenic environmental stressors, we project salt intrusion within the Mekong mega-Delta, in the next three decades. We assess the relative effects of various drivers and show that anthropogenic forces such as groundwater extraction-induced subsidence and riverbed level incisions due to sediment starvation can increase the salinity-affected areas by 10–27% compared to the present-day situation, while future sea level rise adds another 6–19% increase. These projections provide crucial input for adaptation policy development in the Mekong Delta and the methodology inspires future systemic studies of environmental changes in other deltas. Human activities, such as groundwater extraction and sediment starvation, are projected to add to climatic factors like sea level rise to exacerbate saline water intrusion into the Mekong Delta, Vietnam over the next 30 years, according to process-based model simulations.


Environmental drivers of SWI
The monsoon-driven seasonal variability of the Mekong River Basin (MRB) is the mainstay of livelihood along the river. In principle, during the wet season (July-October), the SW monsoon floods the basin and fills the Tonle Sap Lake (TSL) in Cambodia (See Fig. 1). During the dry season (Dec-May), the strong NE monsoon wind increases water levels along the VMD coast (see Fig. 1b), which together with low discharge and higher tidal amplitudes 58 drives salinity intrusion up the delta. Total dry season freshwater discharge is nearly 5-10% of the yearly discharge in Kratie, Cambodia (see Fig. 1a). Historically, during the dry season, the pristine MRB provided two sources of freshwater to the VMD: (1) the Mekong River and (2) the TSL that was filled in the wet season and acted as a flood retention area. While the MR dry season discharge has increased over the past three decades 13 (driven by CC and hydropower operations), due to trimming of the dry season flood pulse by the dams, nowadays, water levels hardly reach the threshold of filling the TSL during the wet season 35,59 and the wet season water levels (Supplementary Fig. 2) follow a declining trend in the past two decades (See methods and supplementary note 1). Therefore, the delta is short of TSL as its second freshwater source in the dry season and experiences a longer period of exposure to SWI. On top of that, as discharge to the delta is not increasing despite the dry season increase of the Mekong River, driven by deeper riverbeds and relative SLR, the delta is much more vulnerable to SWI 60 . Because of this combination, actual SWI in 2016 and 2020 already surpassed 61 the 2050 projections, developed in 2015 26 .
The CC-related drivers are mainly reflected in SLR and upstream discharge anomalies 62 . Fig. 2a shows the projected SLR 63 along the VMD coast. Fig. 2b shows trends of cumulative seasonal upstream discharge variation, based on the Integrated Catchment Model (INCA) 55,64,65 , forced with downscaling precipitation and temperature for two representative concentration pathways (RCP) of moderate 1.5-2°C (RCP 4.5) and extreme 3-4°C (RCP 8.5) global temperature rise 55 . These trends were validated by a thorough comparison with trends simulated by PCR-GLOBWB 66 under different climate forcing scenarios used in the Aquaduct Flood Project 66,67 (see "Methods" for further detail). To provide a realistic projection of discharge, we chose the time series of an average (non-drought) year (in this case 2015) as baseline, and uniformly increased the time series. This increase was calculated for the cumulative discharge of the dry period to match the projected temporal trends. In this manner, we account for long-term upstream discharge variation, but also maintain a typical intra-seasonal discharge variability.
The anthropogenic drivers are categorized under extractioninduced land subsidence and riverbed level changes due to sediment starvation. Fig. 2c shows two opposing scenarios of groundwater extraction 57 , leading to projected land subsidence 57 due to aquifer-system compaction. With rising awareness on the consequences of groundwater extraction, a scenario that reflects 5% annual reduction to a stable 50% of the assumed total extraction volume in 2018 (M2) is considered against a business-as-usual scenario (B2) with 4% (of 2018) annual increase in extraction, similar to the highest rates witnessed in the last 25 years 57 . These subsidence rates are translated to local bed level changes of canals and rivers. Furthermore, over the last two decades, the estuarine system has deepened on average 2-3 m 13,68 , equivalent to losing 2-3 G m 3 (150-200 M m 3 yr −1 ) due to sediment starvation in response to upstream sediment trapping and downstream transboundary sand mining 13,34 . With foreseen 1 G m 3 sand demand until 2040 69 and continuous erosion due to upstream sediment trapping, we consider two trends until 2040 (Fig. 2d). Identical erosion rates as the past 20 years under RB3 (business-as-usual) and significantly lower (one third) erosion rate until 2040 (RB1). RB1 is motivated by rising awareness, shortage of erodible material and potential policy changes among the LMB countries. Given the wide range of possible combinations, we grouped the moderate and extreme (business as usual) scenarios separately to provide a range of potential future projections (also see "Methods").
Results: future of saline water intrusion A recently developed, state-of-the-art, 3D numerical model 48,58 of the VMD (Delft3D-FM) (Fig. 1), stretching from Kratie to 70 km offshore, simulates salinity in the main estuarine channels and the network of primary and secondary irrigation canals over the entire delta. Fig. 1 shows the wide range of possibilities for the spatial distribution of maximum (P100) future SWI in the VMD. The impacts of individual and combined drivers are summarized in Fig. 3 that shows spatial variation of 50th percentile (P50) net salinity increase, relative to a present normal year (2015) and the  57 and riverbed level incision due to sediment starvation (RB1 & RB3) laid over the digital elevation map of the VMD; 93 RCP4.5 + M2 + RB1 (moderate scenarios all drivers) defines a scenario with climate change-driven discharge variation and sea level rise combined with M2 land subsidence and RB1 riverbed level incision (likewise for scenarios all drivers); the surface water numerical model domain and the provincial map of the VMD (coord. system WGS84-UTM 48 N); Monthly variation of cumulative discharge in Kratie with mean and the envelope in gray (a) and monthly averaged water level in Binh Dai with the envelope in gray (b); Estuarine branches (c).  Fig. 3a-f), can increase salinity on average by 2-2.5 ppt by 2040 in all coastal provinces and add~6% additional salinity-affected areas. However, within the estuarine networks, the effect of SLR is balanced by upstream discharge increase ( Fig. 3a-f). Similarly, additional land subsidence ( Fig. 3g-l), especially under extreme scenarios, mainly exacerbates SWI in the coastal provinces (additional average 0.5-1 ppt by 2040), and adds 2-6% additional areas affected by SWI with limited effect on SWI within the estuarine channels. Riverbed level changes drastically increase SWI in the estuaries (5-30 km further inland) and consequently in large parts of the delta further inland. Even moderate sediment starvation scenarios (RB1) portray a future that under normal circumstances salinity intrudes up to 100 km inland. Bed level changes, depending on sediment supply and sand mining practices (within and beyond the delta) can increase the total affected areas by additional 10-25% by 2050.
The Ca Mau Peninsula shows most vulnerability to SLR and land subsidence. However, similar to the observed trends in the past two decades 13 , the all drivers scenarios show stronger response to the drivers in the Tien River estuarine system (Northern branch) than in the Hau River (Southern branch). This means that not just the coastal provinces, but also more inland provinces, such as Vinh Long and Dong Thap that presently do not experience SWI, will potentially have limited/no direct access to fresh river water during the dry season, and potentially need diverted freshwater from even further upstream or resorting to different types of storages. Under those scenarios, in an extreme drought event, salinity can intrude up tõ 150 km inland (see Fig. 1 and Supplementary Fig. 5-6). Note that the difference in SWI between Tien and Hau Rivers is mainly driven by the fact that in contrast to the Tien river, we expect limited/no change in the riverbed levels of the Hau River, downstream of Can Tho 70 . Sensitivity analysis demonstrates that if the bed levels do change, this has a significant effect on SWI, e.g., for the Can Tho province (see Supplementary Fig. 5-6), but the overall conclusions remain unchanged. The all drivers scenarios demonstrate the shear impact of bed level changes on salt intrusion in areas where saline water never occurred before, raising alarm for more upstream provinces.
The exact effect of CC on SWI is in the interaction of upstream discharge anomalies and downstream SLR. Fig. 4 summarizes the expected variations in SWI-affected areas in time. Under scenarios that only account for CC, the increasing upstream discharge trends until 2060 (see "Methods") results in a slight decline of SWI between 2040 and 2050. However, this cannot be extrapolated to the second of the century, as the projections diverge post-2060. The effects of land subsidence are more local, and especially impact already saline areas with a smaller impact on the aggregate of SWI-affected areas. Riverbed level changes show a substantial impact on spatial distribution of SWI. Note that flattening all drivers curves between 2040-2050 is mainly due to the assumed limitation of riverbed level changes past 2040. All these estimates assume that SLR follows the 5 th IPCC assessment report (AR5) 63 . Considering uncertainties in SLR projections, the sensitivity analysis (Supplementary note 5) shows that an extreme SLR of 60 cm by 2050 in response to Antarctic ice-sheet melting 71 can add an additional 9-13% to the aggregate SWI-impacted areas, with dramatic implications for land use and adaptation strategies (See Fig. 4 and the Supplementary Fig. 9-10).

Discussion
Climate adaptation studies of a mega-delta require in-depth integration of knowledge from natural systemic response to upstream and downstream climatic and anthropogenic developments. Salt intrusion, a key determinant of land use, has long been considered as a byproduct of climate change and specifically global SLR, neglecting the effects of local land subsidence and riverbed erosion. Here, for the first time, we distinguish the effects of main "drivers of exposure and vulnerability" 30 in a mega-delta (for the case of Mekong Delta) and identify that while relative SLR impacts salinity in parts of the delta, riverbed level changes are the greatest threat to other parts of the delta. It is the combined effect of stressors that draws a dramatic image of the future. Unanimous assumption and hence, the conclusion of the previous studies 26,27,49,72 was that CC and SLR are the main drivers of SWI. In addition, they also relied on 1D or 2DH numerical models, excluding vertical variable density dynamics, and limited their domain coverage (see Supplementary note 2 and  Supplementary Table 1). Here, we included the entire delta and extended the model upstream to where tide is diminished and downstream to where salinity is considered constant to minimize errors in flow division within the multi-channel estuarine system 58 and appropriate baroclinicity in river-ocean interaction 73 . Leveraging contemporary numerical techniques and computation power, the application of a high-resolution 3D model captures the 3D dynamics of SWI and the spring-neap estuarine variability as a key salt transport mechanism 48 . Our results disclose that global climate change, given the uncertainties in Antarctic ice-sheet melting can lead to 5-19% increase in SWIaffected areas by 2050. Subsidence can increase this by 5-6%, with more local impact on salinity in the coastal provinces. However, the largest increase in SWI is driven by riverbed level changes with even moderate scenarios resulting in tens of kilometers more SWI and potentially 10-30% increase in SWI-affected areas by 2050.
These findings in projected SWI have substantial implications for the Mekong and other deltas around the world (e.g., Ganges Brahmaputra, Ayeyarwady, Pearl, etc.). While CC threatens the world's deltas, anthropogenic drivers, largely reflected in sediment starvation and hydrological regime shifts, seem to outpace CC in the first half of the 21 st century; therefore, some of the existing trends can be partially mitigated by domestic measures such as strict groundwater and sediment budget management.
Our integrated systems approach, incorporating regional, local, natural, and human-induced drivers highlights the sensitivity to often uncertain socio-economical or geopolitical drivers, and underscores the invaluable role of preserving delta elevation and riverbed levels in deltaic systems.
In the case of the Mekong Delta, the trends are clear. Upstream entrapment leads to bed, bank, and coastal erosion, and alters the hydrology cycle, disrupting the crucial retention functionality of the TSL. Due to deeper channels, the delta itself is far more vulnerable to change (e.g., drought); sand mining further deepens the rivers and groundwater extraction lowers the rest, all inviting salt further inland. Reducing groundwater extraction and sand mining can save up to 600 10 3 Ha of land from saline water. However, as encroaching SWI strains freshwater supply, shifting from non-renewable groundwater to surface water demands a delicate strategy. Massive socio-economic implications of land use transformation at such scales, may disproportionately affect the vulnerable communities in and beyond the delta. This calls for enforcing policies that at least maintain the status-quo, as given the sediment deprivation of the Mekong and other big river deltas 11 , returning to the better past seems inconceivable. While swift enforcement of effective policies may portray a possibly sustainable future for the deltas, inertia in translating new policy to effective mitigation measures can lead to further irreversible degradation of the vulnerable freshwater systems and their associated socio-economic structure.

Future work
The present work is a multi-disciplinary study, combining efforts of hydrologists, coastal oceanographers, and geologists, and results of the seven-year international "Rise and Fall" research project. It forms the foundation for potential follow-up work and future improvements, in particular on the following aspects: (1) its time span, (2) processes, (3) data, and (4) adaptation and mitigation strategies. Our current study is limited to the first half of the 21st century (until 2050). Beyond 2050, the global and regional climate models, and thus climate scenarios diverge, potentially leading to a wider range of pathways for SWI. In terms of processes, the present work does not include (a) evaporation/ transpiration variations driven by climate change in the 3D deltawide model, (b) variations in monsoon winds in the next decades, (c) tidal dynamic variations in response to SLR (more important beyond 2050), and (d) exact morphological response of the system to sediment starvation. The latter is an extremely complicated affair as it must combine a sophisticated monitoring campaign (to collect bed composition in depth), derive an accurate account of sediment budget in the delta, including sand mining within and beyond the delta, to be able to carry out a reliable modeling exercise. Regarding data, the first step would be to update climate forcing with CIMP6 data for upstream discharge simulations and incorporate accelerated SLR based on the sixth IPCC Assessment Report (AR6). The present manuscript is a first attempt to assess potential consequences of reducing water demand, and different scenarios of groundwater extraction-induced subsidence. Future work on exploring these morphological interactions and feedbacks of human water use and management of subsidence, sediment dynamics and SWI holds promise to deepen our system understanding and potentially work towards integrated solutions. The present framework can be employed to test solutions and assess/evaluate the consequences of various adaptation measures in time, aiding the development of adaptation pathways for the VMD.

Methods
Modeling and data. We applied a state-of-the-art high-resolution 3D surface water numerical model, extended from 300 km upstream to 70 km offshore, coupled with a 3D geological model of the VMD to project salt intrusion. We developed boundary conditions from various global climate models and is discussed in next paragraph. The estuarine channels were modeled in 3D, with horizontal grid resolution of 200-250 m and vertical z-layer resolution of 0.5 m. The upstream channels and the irrigation system of primary and secondary channels were modeled in 2DV in similar horizontal and vertical grid resolution. The 3D model domain ensures correct advection of solid matters as opposed to the commonly-used 1D modeling approach of estuarine channels that translates non-linear salt transportation mechanisms 48,74-79 such as stokes transport, tidal pumping, tidal straining, and gravitational estuarine circulation to a single dispersion coefficient 80 . The model bathymetry of the primary and secondary channels was provided by SIWRP from 2008, but the river and estuarine bathymetries were updated with the 2018 collected survey data 68 , which guarantees an accurate starting point for future projections. The gauge data within the VMD (water level, discharge, and salinity) was collected by Southern Regional Hydro-Meteorological Centre (SRHMC) and received from Southern Institute for Water Resources Planning (SIWRP). The existing and planned sluice gates were incorporated in the model (see Fig. 5).
The 3D surface water model of the VMD 48,58 was forced by daily averaged discharge from upstream (Kratie, Cambodia) and water level from downstream offshore boundary with assumed constant salinity of 35 PSU. The water level boundary condition consisted of 13 leading tidal constituents (M 2 , S 2 , N 2 , K 2 , K 1 , O 1 , P 1 , Q 1 , MF, MM, M 4 , MS 4 , MN 4 ), the subtidal ocean surge and the average sea level (0 m at present and increasing in time to replicate SLR). The tidal components assumed to be constant throughout future scenarios, but we tested the model sensitivity to typical range of variations in our sensitivity analysis (supplementary Fig. 7, 8). The subtidal ocean surge, similar to previous work by the authors, was superimposed to the offshore boundary with a phase lag from the measurements at the estuary mouth 58 . We modeled existing water demand 58 as a sink and considered the existing water demand figures to not change in the future. However, given the envisioned reduction of water demand in the future, we tested model sensitivity to 25% lower extraction. We assumed that evaporation increases from 4 to 7 mm per day from January to April 58 (as present day) with no precipitation during the dry season.
The numerical domain consists of 3D (the sea and the estuarine channels) and 2DV (width-integrated but depth-varying in the system of primary and secondary channels) sub-domains. Horizontal resolution is~200 m along-channels and 50-200 m across-channels (depending on channel width). In vertical, the numerical domain follows a so-called z-layer scheme with 0.5 m resolution (30 layers in 15 m depth). Vertical resolution gradually increases below 15 m depth in the deeper continental shelf. Following 58 , the re-calibrated spatially varying Manning bed friction coefficient is lower along the downstream estuarine distributary channels (0.02 s m −1/3 ) and increases to 0.04 s m −1/3 along the upstream rivers. This not only results in a good comparison of water level and discharge against the observations, but also follows the observed riverbed composition, with muddy seabed near the sea and increasing sand content when traveling upstream 81 that is associated with higher bed roughness. The numerical model was extensively calibrated 48 for the storm event of the dry season of the year 2016, and in this study, we validated the model for the dry season of the year 2015.
Seasonal model simulations were carried out between January 1 and April 30. The first 30 days of the simulations were considered as spin-up time and the postprocessing was done over the model data of February, March, and April. Fig. 6 presents an overview of the model performance for the year 2015. The model shows excellent quantitative performance for water levels and discharges across multiple stations. Given the complexity of comparing measured and modeled salinity 48 , the model performs reasonably well in capturing salinity magnitude, trends, variation, and especially the spatial extent (e.g., see Fig. 6d, k & l). Previous publication by the authors expands extensively on the background, development, and performance of the applied model 48 . The delta-wide spatial variation of modeled salinity for the existing condition is presented in Fig. 5. Note that also in the existing conditions, for sake of consistency with the future conditions, the planned (under construction) salinity sluice gates are already included. Any variation in the projections is compared to this simulation, representing an average year salinity distribution in the dry season (2015).
Environmental drivers. Effects due to CC are incorporated in upstream discharge projections 55 and SLR 63 for RCP 4.5 and 8.5 scenarios. There are uncertainties with respect to both parameters; however, the consensus around increasing sea levels infers strong evidence that SLR is the main driving force of SWI [26][27][28][29] . On the other hand, while there is also a consensus in increasing dry season discharge due to climate change and hydropower operations [51][52][53][54][55] , the effect of this opposing factor on SWI was undetermined. As this is a positive driver, in order not to underestimate SWI, while applying published discharge projection trends 55 , we validated  83 and withdrawal for irrigation, livestock, industry, and households. If no projected data were available, the model fell back to the latest available data point in time. After accounting for water use and withdrawal, excess precipitation is routed using the kinematic wave simplification of the Saint-Venant equations. For further information, consult the model description paper 66 . Using a different model and data combinations, any interference between the INCA data 55 and the validation data set can be avoided, bolstering the applied trend. Note that the INCA model is a Regional Climate Model (RCM) and the PCR-GLOBWB is a Global Climate Model (GCM); therefore, there can be differences (offset) in their computation of Mekong River discharge. Nevertheless, the trends should not demonstrate bigger differences, as they do agree in this study.
Since running high-resolution 3D models continuously for multiple decades is neither smart nor feasible, instead of the absolute projected discharge values, we initially derived increasing trends of projected dry season cumulative discharge. Then, by gradually increasing the 2015 dry season discharge time series, so that cumulative discharge matches the projected trends we translated a present-day normal year (2015) to the next decades (see Fig. 7 for cumulative discharge trends). Following this method, while committing to natural existing conditions, we refrain from exaggerating future SWI. This may disregard intra-seasonal variability of discharge in the coming decades due to climate change. However, given the fact that flow condition will be highly regulated by hydropower operations, we expect little to no impact on the projections. Furthermore, in projecting upstream freshwater discharge, we have excluded the freshwater supply of Tonle Sap Lake during the dry season. Historically, Tonle Sap Lake (TSL) provided 25-30% of dry season freshwater supply of the VMD, and generally fed the delta between November and February 58,84,85 . After February, the sole freshwater supply to the delta was always the Mekong River. In recent years, due to hydropower operations and smoothing of the flood pulse, the duration of the dry-season freshwater supply from TSL has decreased. During the 2019-20 extreme drought of the Mekong, already in November, the TSL could not supply freshwater to the VMD. Therefore, as we focus on the dry season regime (Jan.-Apr.), we assume no freshwater supply from the TSL during the dry season for the projections.
The two major anthropogenic changes incorporated in this study are (1) land subsidence due to groundwater extraction 42,43,57 and (2) riverbed level incision 13,44,46 due to sediment starvation, driven by upstream impoundments and sand mining. Two opposing scenarios of extraction-induced land subsidence were selected 57 , one business-as-usual scenario in which groundwater extraction continuous to increase (B2) and on mitigation scenario in which groundwater extraction is decreased ( Table 1). The bed levels of the estuarine and irrigation system of primary and secondary channels are adjusted based on linear interpolation among recently developed cumulative subsidence projection maps 57 . The two opposing scenarios of riverbed level incisions (Table 1) are chosen based on the recent estimates 13 of average riverbed level changes over the past 20 years. Riverbed topography of the VMD is not only a function of tidal and flow variations in different time scales, but also a function of upstream hydropower development and sand mining. Based on the past trends, we developed two scenarios: RB1) Until 2040, riverbeds incise one third of the last twenty years and RB3) until 2040, riverbeds incise similar to the last two decades. If policy development leads to more favorable sediment budget in the next three decades, it can be deduced in other projection maps. Note that we assume for both scenarios that riverbed level incisions do not progress beyond 2040. This is motivated by raising awareness in time, and limitations in erodible/available material at the riverbeds.
Several other assumptions were made regarding the environmental drivers. In projecting future elevation change of the delta system due to extraction-induced subsidence, we assumed the digital elevation model of the delta 28,45 to be the true elevation relative to local mean sea level (2018 for riverbed levels and 2008 for the irrigation system). Additional elevation loss as a result of natural compaction of Holocene sediments 86 was assumed to be counterbalanced by deposition of new sediments and therefore omitted. Based on significant scientific evidence 13,70 , we assumed that riverbed levels in the Hau River, downstream of Can Tho will remain stable. Because of its potential significant impact, sensitivity of the system given this assumption is addressed in the supplementary information. We did not study any potential effects of changes in monsoon wind patterns and perhaps variations in average water levels along the Mekong coast during the dry season. Our early assessment did not flag significant trends in subtidal water levels along the Mekong coast over the past two decades other than those associated with SLR 13 . Table 1 summarizes all drivers considered in this study as well as an examined sensitivity analysis (see supplementary information). Furthermore, the simulation IDs are provided as a guide to the repository where the projections can be downloaded 87 .

Time series analysis
Salinity percentile maps. PXX percentile defines salinity concentration that was not exceeded during XX% of the simulation time. To derive spatial variation of salinity over the VMD, the 3D map-files, containing hourly time series of salinity at every vertical and horizontal grid cell over the entire numerical domain was analyzed. Every simulation contained~500 GB of data. The dry-season time series of computed depth-averaged salinity at every horizontal grid cell was used to calculate various percentiles (e.g., P50 and P100/maximum salinity).
Cumulative discharge. The Tonle Sap Lake is filled during the wet season, when the Mekong River's flow partially diverts through the Tonle Sap River to the lake.   During the dry season, the Tonle Sap River flow reverses and drains the lake to the Mekong River and contributes to the total freshwater supply of the VMD (see Fig. 8). To analyze cumulative flow distribution, within the VMD, we first de-tided the discharge signal, by applying the low-pass Godin-filter 88 . Similar to previous work by the lead authors 13,58 , subtidal discharge was then integrated over the dry season. At Kratie, daily non-tidal discharge was simply integrated over a day. Due to lack of data from Cambodia, in order to derive TSL contribution to the VMD, the dry seasons cumulative discharge of Kratie was subtracted from the sum of cumulative discharge in Chau Doc and Tan Chau.
Trend analysis. nonparametric Mann-Kendall test and a Theil-Sen estimator 89,90 were used to fit linear trends to the presented time series. The estimator takes the median of all pairwise slopes of values in time. The method is robust, i.e., it is not sensitive to leveraging or influence due to potential outliers and does not require the assumed normality of the residuals that are associated with ordinary least squares regression modeling. It is suitable for discharge time series that are influenced by various processes 90 (e.g., dam regulations).

Data availability
Delft3D-FM is an open source numerical model https://oss.deltares.nl/web/delft3dfm. The underlying gauge data (observed water level, discharge, and salinity) provided by the SIWRP, following the organizational policy, can be provided upon request for noncommercial use from a limited access repository 91 (https://doi.org/10.5281/ zenodo.4771261). The 2015 wind data originally developed by NCEP/NOAA, can be downloaded from the DHI repository under https://www.metocean-on-demand.com/ #/main. Projections of upstream discharge from the five GCM's and downstream sea level rise can be openly accessed in an online repository 92 (doi: 10.5281/zenodo.4771240). For the spatio-temporal variations of land subsidence 57 we refer to the original publications and its lead author philip.minderhoud@wur.nl. The geometry data of the network of primary and secondary channels, following the internal policy of the institute, can only be provided in direct communication with the SIWRP (siwrp@siwrp.org.vn). The 2018 contemporary bathymetry data can be downloaded from https://hydra.hull.ac. uk/resources/hull:17952. The salinity projection maps in the next three decades and the associated sensitivity analyses are available to be downloaded from an open-access repository 87 (https://doi.org/10.5281/zenodo.4772967).

Code availability
The pre-/post-processing codes used in this study are all written in open-source Python 2.7. Access to specific parts of the codes will be granted on request to the corresponding author.
Received: 1 March 2021; Accepted: 10 June 2021;  Summary of the drivers of environmental change, list of simulations including those for sensitivity analysis.