Modelling Holocene analogues of coastal plain estuaries reveals the magnitude of sea-level threat

Hydrodynamic modelling of Australia’s lower Murray River demonstrates the response of a large coastal plain estuary to the mid-Holocene (7,000–6,000 yr BP) sea-level highstand. The approximately two metre higher-than-present sea level during the highstand forced the estuarine limit upstream generating an extensive central basin environment extending more than 200 kilometres from the river mouth (143 kilometres upstream of the modern tidal limit). The geomorphic history of the region does not conform to conventional estuarine facies models as, for much of the Holocene, the lower Murray River acted as a landward, gorge-confined extension of the Murray estuary. The incredibly low relief of this coastal plain system drove significant saline incursion and limited current velocities across the estuary facilitating deposition of a laminated silt-clay sequence which our results suggest may be regionally extensive. Variations to discharge, barrier morphology, or the estuary’s bathymetry result in minimal change to the estuarine palaeo-environment. The shift to the present-day fresher water distribution in the Murray estuary requires a fall in sea level to present-day conditions. The dominance of sea level as the controlling factor on this estuarine palaeo-environment highlights the significant potential impact of climate change induced sea-level rise to coastal plain estuaries.


Supplementary methods: Model calibration
Following a significant period of prolonged drought, BMT WBM were commissioned by the Murray-Darling Basin Authority (MDBA) to perform a feasibility study assessing the adoption of a virtual weir at Wellington (rkm 78) to provide an adequate fresh water supply for the region. This involved the establishment, calibration and validation of a hydrodynamic model encompassing the LMR and Lower Lakes region downstream of Blanchetown (rkm 282, Lock 1). A validation dataset was collected, and the model calibrated for hydrodynamics -including water level, wind and river fluxes -and salinity. Model calibration adequately resolved short term (i.e. a single extreme saline intrusion event) and long term (i.e. entire 17 month dataset) trends.
Owing to the long-term temporal scale of modelling estuarine response to sea level change over the course of the mid-to late-Holocene, simplifications need to be applied based on best estimate assumptions to guide parameters for sensitivity testing. The model set up adopted in this study seeks to apply 'appropriate complexity' balancing a reductionist approach to input data based on geological correlation, to produce outputs which are computationally efficient yet meaningful (1, 2). Model manipulation and model scenarios involved a deviation from present day morphology, flow and flow obstructions, sea level and ocean outlet, resulting in models which inherently could not be calibrated against the present day. Instead, results were compared to the Holocene stratigraphic record. Water heights were compared to documented evidence of notches and wave-cut cliffs along the former shoreline of Lake Alexandrina (3), and inundation extents correlated with the Malcolm soil combination and sediments of the Cooke Plains Embayment (4,5).
Wave data was excluded from our model as the primary influence of waves within a wave-dominated estuary is as a driver of morphological change through the formation of a barrier complex at the estuary mouth and this model does not incorporate a sediment transport or morphology component (6). Furthermore, although this estuary is wave-dominated at its entrance, owing to the immense scale of this system, areas subject to significant wave energy present a negligible component of the overall model domain (Fig. S5).
A Manning's coefficient of 0.025 was adopted across the model domain for this study. Applying a varying Manning's coefficient was not actually implementable without a robust understanding of surficial sediments at the Holocene highstand over the entire 282 rkm of the model domain. Given this impracticality, applying a global Manning's coefficient was deemed sensible and, although not an accurate representation of reality, this method nonetheless provides the means for direct comparison between results. The value of 0.025 was selected as it lies within the bounds of appropriate Manning's coefficients given the likely Holocene palaeo-environmental conditions (7). Further, this value was deemed appropriate as sensitivity testing conducted during calibration of the base model provided by BMT WBM revealed that results did not vary significantly given changes in the Manning's coefficient, but the model was best resolved when adopting values between 0.015 and 0.02538 (8). Reference 8 is a government funded study. For access, please contact the Murray Darling Basin Authority (MDBA) through: engagement@mdba.gov.au.

Supplementary methods: Morphology
To best assess the influence of morphology on the hydrodynamics of the system at the Holocene highstand, three surfaces were created. A pre-regulation surface (S up) provided a modern-day end member, the depth to the Monoman -Coonambidgal Formation transition provided a late-Pleistocene -early-Holocene end member (Slow), with the third surface a best estimate of highstand bathymetry and topography (Smid). The Slow surface is certainly deeper than at highstand, with the average depth of the body of Lake Alexandrina (between Point Sturt/Point McLeay and Pomanda Embayment) approximately -43 m AHD (Australian Height Datum; approximately mean sea level). By comparison the Sup surface has an average of approximately -3 m AHD over the same area, with the Smid best estimate highstand surface at approximately -8 m AHD. Within the LMR, the average of three surfaces is more closely constrained, varying from approximately -15 m to -8 m.
To resolve the Sup pre-regulation surface, the lock, barrages (and associated sediment sills), man-made levies and modern flood tide delta were removed from modern day DEMs (Fig. S6a). To resolve the Slow Pleistocene -Holocene surface, depths were interpreted from over 100 sediment cores, as well as interpretation of data and maps by Barnett (9) and Von der Borch and Altmann (10) (Fig. S6b; Table S3). The location and depth of the palaeo-channel within Lake Alexandrina was based on an interpretation of the work of Barnett (9) and geological maps (Fig. S6b). The Smid best-estimate highstand surface has the greatest uncertainty as it was resolved by subtracting regional sedimentation rates from the pre-regulation surface (9,11) and dated sediment cores (n = 18; Fig. S6c). Within the LMR and thalweg seaward, a sedimentation rate of 0.69 mm/y was adopted (9,11). All other elements were adjusted with a sedimentation rate of 0.16 mm/y (9,11).
The spatial interval of the three bathymetric surfaces created is identical to that of the modern-day input dataset as values were adjusted at each individual cell. The cell/element size varies considerably across the model domain with an average of approximately 50 m in cell side length within the LMR, modern flood tide delta (seaward of Point Sturt/Point McLeay, rkm 40) and Lower Lake fringes, and palaeo-Murray thalweg, with a considerably larger cell/element size within the main body of the Lower Lakes.
The chain-of-islands evolution of Sir Richard and Younghusband Peninsulas (3,4,12,13) was the guiding premise for the series of barrier morphologies presented in this study. Bourman and Murray-Wallace (13) and de Mooy (4) give detailed descriptions and maps of hypothesised former outlets of the LMR to the ocean. These maps were georeferenced and digitised and assessed relative to the modern-day topography of the barrier system. Furthermore, the presence or absence of Aboriginal middens within the Holocene barrier (3,14,15,16) were also mapped and assessed relative to their radiocarbon ages. These data were combined and analysed to produce a series of best estimates of the morphology of the mouth of the LMR as it evolved throughout the Holocene. As this study imposes a static barrier configuration for each scenario, an assessment of the dynamic response of the barrier to events on short temporal scales, such as tides or storms, was beyond the scope of the study.

Supplementary methods: Comparison of 2D and 3D simulations
Estuarine stratification can cause a salt wedge at depth which cannot be resolved by a 2D simulation. The presence of a salt wedge has the potential to alter inferences drawn from estuarine zonation and the likelihood of deposition of a laminated sequence whereby flocculation may be assisted by salinity. The scale of the study area, with the model domain spanning some 282 river kilometres, precluded the use of a 3D model setup without justification as the computational power to run such a simulation is ten times that of its 2D counterpart. Therefore, a representative subset of models were run in 3D to assess the suitability of adopting 2D models for this study. This representative subset of models allowed for a comparative assessment of 2D and 3D results for all bathymetric surfaces (S lowWL2DavBmod, SmidWL2DavBmod and SupWL2DavBmod), both sea level scenarios (SmidWL2DavBmod and SmidWL0DavBmod), all discharge scenarios (SmidWL2D-Bmod, SmidWL2DavBmod and SmidWL2D+Bmod) and all barrier morphologies (SmidWL2DavB0, SmidWL2DavB+, SmidWL2DavB++ and SmidWL2DavBmod).
3D simulations retained the same model setup with the inclusion of a parametric vertical mixing model with a second order vertical solution and density coupled salinity. A comparative analysis of results suggests that 2D simulations are sufficiently representative of 3D simulations (Fig. S7). Minor differences in maximum salinity reached across the region does not alter the designation of the palaeo-environment with 2D simulations providing a conservative approximation of 3D results (Fig. S7 a-b). Crucially, the brackish limit (equivalent to 1 psu) differed by a maximum of 1 rkm between 2D and 3D depth averaged results, with the exception of the pre-regulation (Sup) surface where the 3D simulation captured a salt wedge that penetrated 6 rkm further upstream (SupWL2DavBmod scenario). A negligible change in the total area conducive to the deposition of a laminated sequence was also observed, with 2D simulations overestimating this area by an average of 4%, demonstrating that 2D approximations of 3D velocity magnitudes are appropriate (Fig. S7 c-d).  (17). Refer to Table S1 for scenario descriptions.   (19,20,21). Refer to Table S1 for scenario descriptions.    Table S3. The location of the palaeo-Murray thalweg (blue) was inferred from data presented in Barnett (9). (c) Creation of the Smid surface was resolved by subtracting regional sedimentation rates from the pre-regulation surface (9,11) and dated sediment cores (red). Within the LMR and palaeo-Murray thalweg seaward (blue), a sedimentation rate of 0.69 mm/y was adopted (9,11). All other elements were adjusted with a sedimentation rate of 0.16 mm/y (9,11). Figure S7: Comparison of SmidWL2DavBmod 2D and 3D key outputs of maximum salinity and velocity magnitude. Maximum salinity reached in a best estimate Holocene highstand (a) 2D simulation is comparable to that of a (b) 3D simulation such that the classification of estuarine zonation remains consistent. 2D simulations provide a conservative approximation of 3D results. Salinity is measured based on the classification scheme of Tooley (17). Maximum velocity magnitude in a best estimate Holocene highstand (c) 2D simulation and a (d) 3D simulation. Areas are shaded red where maximum velocity > 0.3 m/s and therefore is not conducive to the deposition of a laminated silt-clay sequence (19,20,21). A 3% change in total area demonstrates the negligible difference in velocity magnitude between 2D and 3D simulations.
Table S1: The scenarios adopted in this study. Codes identify bathymetric surface (Sup = premodification condition, Smid = highstand best estimate condition, Slow = Pleistocene-Holocene boundary condition), discharge (D-= drought, Dav = pre-regulation average, D+ = pre-regulation average with a flood event), and barrier morphology (B0 = no barrier, B+ and B++ = phases of chain-of-islands evolution (3,4,12,13), and Bmod = present day). A code and category has been assigned to each scenario to facilitate interpretation. Scenario categories are grouped based on sea-level and bathymetric surface. Table S2: Description of stratigraphic formations and soil combinations within the study area. Table S3: Sedimentological data used to inform the creation of the Slow surface. Core logs (9,23,24) and cone penetrometer tests (25) were analysed for the transition from the Monoman to Coonambidgal Formation, which is interpreted to represent the approximate Pleistocene to Holocene boundary (Slow). Depths were adjusted relative to present day bathymetry and topography to give elevation in meters AHD.