The Fitzroy Estuary (FE), located at ~23.5°S (see Figure 1 ), extends from a barrage at Rockhampton, ~70km from the coast and winds through low-lying agricultural catchments, finally discharges into Keppel Bay (KB) via southern deep channel. The FE-KB system in this study includes the upper FE (Barrage to Thompson Point), lower FE (Thompson Point to Sea Hill point, including southern tidal creeks) and KB (north of Sea Hill). Upstream of the barrage, the Fitzroy River (annual range of discharge 0 – 14,000 m³/s) drains a grazing dominated catchment area of about 142,000km². Due to the catchment’s large size, there is a wide range of soil types, land cover characteristics, and rainfall-runoff dynamics, leading to climatic, hydrological, and geological differences between sub-catchments (Dougall et al., 2014; Huggins et al., 2017). Numerous studies have reported high level of sediment and nutrient loads from the catchment to the reef since European settlement (McCulloch et al., 2003; McKergow et al., 2005; Brodie et al., 2011). The current sediment load from the FE-KB catchment is estimated to have tripled due to increased land clearing and changes in land use since pre-European settlement (Kroon et al., 2012). The depositional history of the FE-KB system indicates that sediment from the river has accumulated in and around the lower FE (Thompson Point to Sea Hill) and the lower KB to form a more deltaic-type environment, characterized by extensive floodplain, salt flats, mangroves, and tidal creeks (Brooke et al., 2006; Ryan et al., 2006; Bostock et al., 2007).

The export of sediment from the FR is mainly controlled by highly episodic summer floods during wet months (November to March) and is redistributed primarily by tidal currents within the estuary and bay during the dry months (April to October) (Webster and Ford, 2010). The monthly-averaged hydrograph of FR ( Figure 2 ) at the Gap (130005A, https://water-monitoring.information.qld.gov.au/) shows high seasonal and interannual variability in river flow due to the impact of El Nino and La Nina cycles. Generally, there are one or two above-mean annual flows, with the occasional large flood event lasting several weeks, and low flows for the remainder of the year. In 2022, the average discharge rate was 43 m³/s, with multiple flow events up to 1500 m³/s ( Figure 2 ), indicating a relatively dry monitoring year. To assess sediment fluxes during a wet year, the model simulated the 2011 flood event, which caused severe flooding in January that lasted for about 18 days. This event set a record for the highest annual rainfall in the FR catchment and peak discharge rate exceeding 13,000 m³/s. The monitored FR catchment had an estimated annual total suspended solids load of nearly 7.0 Mt. y^-1 in 2010-2011, with the fifth-highest flood peak recorded at 9.2m in Rockhampton.

The FE-KB system is tide-dominated with tidal range reaching 5m during spring tides and current speeds may exceed 1 ms^-1 near the FE mouth. The system is well-mixed during the dry season, while during high river flow freshwater effectively flush the estuary and lower salinity in KB. It takes approximately 6-8 months for the salinity to increase back to marine levels (Herzfeld et al., 2005). Prior studies have shown a wide range of variations in salinity from the estuary head to the lower FE, which can range from being entirely fresh during river discharge pulses to being hypersaline throughout dry years (Webster and Ford, 2010). Port Alma in the lower FE was identified as the zone of maximum resuspension with a higher percentage of mud in surficial bottom sediments (Ryan et al., 2006; Webster and Ford, 2010).

Long term surface turbidity samples from the past 30 years (^© State of Queensland (Department of Environment and Science) 2019) demonstrated the existence of a Turbidity Maximum Zone (TMZ) between stn.35 and stn.40 (refer to Figure 1 for stns labelled as distance from Rockhampton Barrage towards downstream). The turbidity levels are generally higher during the wet season, with greater mean and range, and lower at the seaward end of the estuary. TMZ is consistently higher during both wet and dry season compared to upstream and downstream regions ( Figure 3 ).

The unstructured finite-volume community ocean model FVCOM, coupled with the UNSW-Sed module (Chen et al., 2003; Wang, 2002) was used to simulate the estuarine circulation and sediment dynamics in the FE-KB. The model domain covered ~50km offshore and ~60km upstream from the estuary mouth towards Rockhampton Barrage. The Rockhampton Barrage acts as a barrier to tidal and salt propagation, indicating the tidal limit of FE. Note that distances in this study are measured as downstream from the Rockhampton Barrage ( Figure 1 ). The grid consisted of 40,930 nodes and 77,730 elements, with variable cell sizes ranging from ~2500m at the open ocean boundary down to 100m in the estuary ( Figure 4 ). The vertical resolution was 20 sigma layers with higher resolution at the surface and bottom layers. The bathymetry data used was the high-resolution (30m) depth model for the GBR obtained from Geoscience Australia (Beaman, 2017).

Open boundary conditions were specified by interpolating the dataset of tidal constituents from the TPX08 global tide model (Egbert and Erofeeva, 2002, https://www.tpxo.net/global/tpxo8-atlas, accessed Sep 2022) to the open boundary nodes. Catchment rainfall-runoff processes were simulated by GBR Dynamic SedNet model framework, which incorporates a paddock-scale pollutant generation model that calculates loads from rivers, side streams and creeks within the GBR catchments (Waters et al., 2014; Ellis, 2018; Waterhouse et al., 2018; McCloskey et al., 2021). The daily river discharge and sediment loads from the FR were added at the upstream boundary with salinity set to zero. The model was spun up for one month before the different model scenario runs with an initial zero velocity field and a uniform salinity field of 27 psu. The bottom drag coefficients is defined as 0.001m for marine bed and 0.0005 m for estuarine bed, time step for vertical diffusion is 0.4s and for horizontal advection is 2.0s.

The UNSW-Sed module was two-way coupled to the hydrodynamic model by allowing SSC to affect the seawater density and bottom drag coefficient (Wang, 2002). The detailed module algorithm can be found in the supplement. The parameters used in the sediment module are summarized in Table 2 . Only one type of cohesive sediment is used in this study, of which the settling velocity is a function of suspended sediment concentration. It includes the flocculation and hindered settling process. Morphodynamics are simulated by calculating changes in bed thickness. The initial bed thickness is set to be 10m uniformly and reduced or increased by sediment erosion or deposition process.

Model scenarios with different river flow conditions combined with tidal forcing in barotropic and baroclinic modes were conducted (summarized in Table 3 ) to investigate the role of discharge magnitude on estuarine sediment transport (trapping and erosion) and the processes controlling the location of the TMZ.

Observational data of current velocity, salinity, and turbidity were grouped into wet and dry season data based on the presence of river flows (> 100 m³/s) within 10 days before data collection. These data were used to reveal the seasonal pattern in hydrodynamics and sediment transport, as well as possible drivers of seasonality in the FE-KB.

Simulated water surface elevation, salinity, SSC, and current velocity were compared with observation data to evaluate model performance, and calibration results are provided in the supplement. The model has the capability to predict the sediment flux along and across estuary within a defined area to provide a closed balance of sediment budget. The sediment flux was determined by integrating the product of water depth, SSC and current velocity at each layer along the cross-estuary or along-estuary transect, from the surface to the bottom. The resulting sediment flux was then separated into two components: mean advection (tidally averaged) and tidal pumping (tidally varying) components (McSweeney et al., 2016; Xiao et al., 2020). To calculate these components, the time series of SSC and principal velocity along -layer were low-pass filtered using a 36hr lanczos filter ( in the subscript is omitted for convenience):

where 〈 〉 t indicates the tidally averaged component and { } t the remaining component; F_adv is the mean-advection component of sediment flux, F_pum is the tidal pumping component of sediment flux, and D is the total water depth. In the along-estuary direction, positive values indicate down-estuary transport, and negative values indicate up-estuary transport. The along-estuary and cross-estuary sediment flux were interpolated along the thalweg line and at the XS1 -XS5 transects at 100m intervals.

To evaluate the contributions by tidal and nontidal processes to the total sediment transport, the contribution (%) of mean advection and tidal pumping was calculated as follows:

Becherer et al. (2016) developed a close balance to decompose F_adv into barotropic and estuarine circulation components. They also used tidal averaging instead of the low-pass filter method, with the M₂ tidal period as the averaging window. In their study, both methods were compared and found to agree well in reproducing the trend of sediment flux component variations, except for an unclosed balance in the previous method due to the omission of cross-correlation terms between mean velocity and fluctuating sediment concentration and vice versa. The estuarine circulation component is further extracted following Becherer et al. (2016):

The components calculated in Becherer et al. (2016) are in tidal-averaged and depth-mean format.

The main shipping channel is maintained at depths of 14m to 18m from the estuary mouth towards Port Alma ( Figure 5A ), resulting in strong tidal influence with high salinity ( Figure 5B ) and intermittent sediment resuspension due to strong tidal currents. Thompson Point (40 km) marks the rapid transition from deep and narrow (~16m deep and ~600m wide) to shallow and wide downstream (less than 10m deep and ~2000m wide; Figure 1 ). Moving upstream from Thompson Point toward the Bend (30 km) tidal motion gradually mixed saline water with freshwater. At the Bend, tidal currents and resuspension effects decreased, as was found in previous studies (Douglas et al., 2005; Webster et al., 2006). The area upstream from Bend to Rockhampton was dominated by the river with close-to-zero salinity and high turbidity throughout the study period. In Nov. 2022, the edge of the salt intrusion during rising tide was observed between Bend and Thompson Point about 32km from the Barrage. A classic two-layer structure of gravitational circulation was observed, with freshwater at the top and dense water at the bottom ( Figure 5C ).

The three mooring stations exhibited similar probabilities of exceeding depth-averaged current speeds, reaching a maximum of 1.5 m/s. The majority of current speeds (over 90%) remain below 1.0 m/s, or with the current speed exceeding 1.0 m/s less than 10% of the time ( Figure 6A ). A1 station stands out with an ebb-dominant asymmetric current, where over 10% of the ebb current exceeds 1.0 m/s during both the wet and dry seasons ( Figures 6D, G ). However, A3 station is located too close to the northern side of the riverbank to represent the main channel flow behaviour, which is not showing the flood-dominant tidal current that previous studies have predicted to induce the up-estuary pumping of sediment (Margvelashvili et al., 2003; Webster et al., 2003; Webster et al., 2006). The across-axis variability in current velocity was confirmed by ADCP transects in Nov 2022 using a SonTek M9 Hydrosurveyor (shown in supplement). In the survey year, A1 and A2 stations were not completely flushed fresh, with 80% of bottom salinity above 30 ( Figure 6B ). Conversely, A3 station experienced strong interactions between river discharge and tidal propagation, with 80% of bottom salinity above 8 and 18 during the wet and dry seasons, respectively ( Figures 6E, H ).

Bottom turbidity during the dry season shows a higher percentage of high values at Port Alma than at the other two stations ( Figure 6F ), though there is a lack of data for Port Alma during the wet season. A1 and A3 stations exhibit similar variation trends, with over 70% of the bottom turbidity during the dry season limited to 100 NTU, while high turbidity events are more common during the wet season ( Figures 6C, F, I ).

Brooke et al. (2006) predicted sediment transport pathways and budgets for fine sediment (<63μm) based on sediment accumulation rates estimated from sediment cores in the FE-KB. The estimated catchment load for annual sediment (4.575 Mt y^-1) used in sediment budget calculations by Brooke et al. (2006) is comparable to that of a wet year (2011) and two orders of magnitude greater than that of a dry year (2022). In Brooke et al. (2006), 55% of the annual average load was estimated to deposit within the estuary, with a large fraction deposited in mangroves (30%) and the lower FE (10%), including southern tidal creeks. The trapping coefficient over the past decades falls within the range reported in Table 3 : 100% for a dry year and 28% for a wet year. Although our calculations indicate that the upper FE is consistently eroding during the model intervals, Brooke et al. (2006) estimated an accumulation of 8.6% within the same compartment. Asymmetric tidal currents were identified as the primary driver for net up-estuary transport (Webster and Ford, 2010), but this phenomenon was only observed within the TMZ and not within the upper FE in our study. Note that the model in this study does not include storage on the floodplain, mangroves and salt marshes, whereas the measured accumulation rates for the FE do include sedimentation in these areas and represent average accumulation over longer timescales. Besides catchment sources, offshore input can also be a crucial source of sediment for a river system, which was not reported in the study by Brooke et al. (2006). Our research has revealed that during low flow conditions, offshore influx can be transported to the system by currents and tides, resulting in a significant contribution to the sediment supply. Following a pulse of sediment export from the system into the offshore zone during a flood event, tidal, wind- and wave-driven currents tend to transport these river sediments back into the estuary where they are deposited in the lower FE and southern tidal creeks.

This study emphasizes the importance of considering the effect of seasonal variability in river flows on sediment accumulation rates when calculating the sediment budget for wet and dry seasons in the FE-KB. Similar modeling practices were conducted by Webster et al. (2006), who analyzed three single-year periods representing low (1993), medium (2003), and high (1999) river flow scenarios in the FE-KB. Both this study and Webster et al. (2006) revealed continuous channel erosion during low-medium flow conditions, however Webster et al. (2006) observed persistent export flux towards the GBR lagoon for all flow conditions. In this study, catchment loads were trapped within the system without export during low flow conditions. The updated sediment budget analysis in this study suggests that during dry seasons with low flow, sediment from catchment sources, channel erosion, and offshore sources tend to remain within the lower FE (64%), southern tidal creeks (23%), and KB (13%) ( Figure 12 ). During wet seasons with high flow, sediment from catchment sources significantly outweighs channel erosion and offshore import. As a result, sediment accumulates in the lower FE (25%), southern tidal creeks (6%), and KB (5%). A considerable amount of sediment (64%) is also exported offshore during this season.

The drawback of calculating sediment budget estimated from sediment cores is its assumption of a long-term mean sediment delivery to the FE-KB system, without considering the variability in climate of sediment inputs. And it’s quite challenging to apply a mean estuarine trapping coefficient estimated from sediment cores to represent the trapping in the estuarine domain of GBR catchments, which is the current practice to determine trapping and erosion from catchment to reef. This study provides a framework for measuring and modeling sediment dynamics along the catchment-to-reef continuum, with the aim of more accurately estimating sediment budgets and estuarine trapping coefficients while taking into account estuarine processes and climate variability. The broader implication is that inputs to marine model from omitting estuarine processes in the catchments makes predicting, and most significantly, comparing the impact of catchment sediment-related water quality impairments tenuous. Livsey et al. (2022) stressed the importance of considering the complex interactions between catchment processes, riverine transport, estuarine dynamics, and reef ecosystem health when studying the catchment-to-reef continuum. They highlighted the differences in estuarine trapping coefficients across various GBR catchments and advocated for more comprehensive and integrated approaches to better understand and manage these ecosystems. By quantifying these interactions, effective strategies can be developed to protect and manage the GBR and other coastal ecosystems.

There are several limitations to this study that should be acknowledged. First, the study period was relatively short, covering only one wet year and one dry year. Future studies could extend the study period to capture interannual variability in sediment dynamics. Second, the study focused only on the effect of tidal forcing and river discharge on the morphological changes of the estuary and did not consider the impact of wind stress and wave-current interactions. Nevertheless, the development of a high-resolution model combined with high-resolution observations in this study advances our understanding of key drivers on sediment dynamics in the FE-KB, which have been previously observed and modelled at relatively coarse spatial and temporal scales.

Results from this study beget several recommendations for further research. First, the use of remote sensing data in combination with field measurements could improve model validation in the spatial domain, as remote sensing can provide snapshots of the entire model domain and even tidal time series in the case of geostationary satellites (Patricio-Valerio et al., 2022). However, in a highly turbid system like FE-KB, optical remote sensing may be limited to surface water and unable to resolve dynamic processes at depth, such as the TMZ observed in this study. Second, revision of long-term monitoring programs to include continuous, high-resolutions measurements at key FE locations, such as near Thompson Point downstream of the deep hole, could provide insights into the long-term trends and variability of sediment dynamics in the estuary. Third, the model is helpful on investigating the effects of future climate change on the sediment dynamics and morphology of the FE. Much of the coastal areas surrounding the FE is composed of low-relief salt flats and floodplains, and high-resolution modelling of sediment dynamics will be critical to predict and manage how this system responds to rising sea levels and more extreme weather in a changing climate. Finally, the model can be used and further developed for decision support, e.g. to inform future dredging and coastal development to promote actions that enhance sediment trapping efficiency within the FE rather than export to the GBR Lagoon.