Hydrodynamic modelling and model sensitivities to bed roughness and bathymetry offset in a micro-tidal estuary

The frequency of hazardous events, such as pollution and floods, are significantly increasing in estuarine and riverine systems due to climate change and human activities. Therefore, monitoring and managing of these water bodies is very important, but providing reliable predictions of material transport in these systems is challenging. Two-dimensional (2D) depth-averaged hydrodynamic models, which are based on the St. Venant equations, are widely used to study the hydrodynamics and particle transport of estuarine and riverine systems (Cea & French 2012; Sandbach et al. 2018; Iglesias et al. 2019a, 2019b). Hydrodynamic models are associated with systematic and random errors, which cause discrepancy between model outputs and observations. Generally, uncertainties in model estimates can result from errors and approximations in measured model inputs (informational uncertainties), model parameters and structure (model uncertainties), and model solution algorithms (numerical errors) (Loucks & van Beek 2017). If modelled variable values differ from the observed values, and nothing randomly occurs, the discrepancy between model outputs and observations can be reduced through a model calibration process to modify model parameters and structure (List et al. 1997). On the other hand, random errors can be optimally reduced by incorporating observations into the hydrodynamic models via a data assimilation (DA) technique (Mazzoleni et al. 2015). Improving hydrodynamic model outputs is important since they serve as the input to simulate the dynamics of particles and sediments in estuaries.

Model calibration can be performed either manually or automatically via adjusting the model parameters. An automatic calibration process is more efficient and reliable. Automated model calibration can be done through various algorithms, such as the Monte Carlo method (Beven & Binley 1992), genetic algorithm (Wang 1991), gradient-based methods (Jiang et al. 2019), and derivative-free algorithm (Garcia et al. 2015). These calibration algorithms are discussed in detail in Efstratiadis & Koutsoyiannis (2010). Bed roughness is usually considered as the main calibration parameter in most hydrodynamic model calibration studies (Piedra-Cueva & Fossati 2007), yet errors in bathymetry also have a significant impact on the estuarine and riverine flow fields (Cea & French 2012; Garcia et al. 2015; Yan et al. 2015). Bathymetry errors can also affect the tidal phase due to the dependence of wave celerity on water depth (Matte et al. 2017). Losch & Wunsch (2003) discussed the importance and practicality of using bathymetry as a control parameter in the ocean circulation model.

Because of the difficulties and costs of conducting bathymetry surveys in some locations as well as inherent errors in some survey approaches, many studies in the literature have proposed that a model can be calibrated by adjusting bathymetry (Cea & French 2012; Zijl et al. 2013; Yan et al. 2015; Wood et al. 2016; Matte et al. 2017; Jiang et al. 2019). It has been shown that by considering errors in bathymetry such as calibration parameters, model performance can be improved noticeably, especially for the estimation of depth-averaged currents for estuaries, which significantly depends on bathymetric data (Wang et al. 2009; Cea & French 2012). Errors in bathymetric data can result from either uncertainties in bathymetry surveys (Cea & French 2012), random errors due to measurement precision, inaccuracies in datum and geographical positioning, or errors that are associated with hydrodynamic modelling (Wang et al. 2009). The latter includes those from the sensitivity of the model to the spatial density of the bathymetric data and the interpolation scheme for mapping bathymetric data into the model mesh nodes. However, there are few studies on the simultaneous incorporation of the uncertainties of bathymetry and bed roughness into the calibration process in riverine (Yan et al. 2015; Wood et al. 2016; Schneider et al. 2017; Jiang et al. 2019) and particularly estuarine systems (Cea & French 2012; Garcia et al. 2015).

During recent decades, the experimental deployment of Lagrangian sensors, such as drifters and floats, has been increasing because of their advantages, including low cost, portability, and more importantly providing larger spatial coverage compared with Eulerian techniques (Strub et al. 2009). In addition, modelling of particles in the Lagrangian frame has been shown to provide more valuable insights to transport of the materials (Suara et al. 2020). Thus, they can provide additional insight into the horizontal motion of particles in estuarine systems. While Lagrangian instruments have been largely used in oceanography (MacMahan et al. 2010; Coelho et al. 2015), there are significantly fewer studies that employed Lagrangian observations in estuarine applications (Tinka et al. 2010; Suara et al. 2015). It is very challenging to incorporate Lagrangian data into hydrodynamic models since Lagrangian observations are measured in a moving coordinate system that are different from model states computed in a fixed frame of reference. However, this type of data offers a new opportunity for DA and validation purposes in estuarine applications (Kuznetsov et al. 2003).

The study region is Eprapah Creek located on the East subtropical coast of Australia and part of the Greater Brisbane area. It has been negatively affected by construction and industrialization, and more importantly, it was significantly polluted in 1998 by illegal discharge of chemical materials (Ferris et al. 2004). In addition, it is a suitable nature's laboratory for researchers due to its small size and low level of recreational activities (Suara et al. 2017). This shallow micro-tidal estuary has been widely studied mostly related to experimental investigation of its turbulence characteristics (Chanson et al. 2005; Trevethan et al. 2007, 2008; Suara et al. 2019), sediment transport (Chanson et al. 2008), water quality (Brown et al. 2004) and the mechanisms responsible for mixing and dispersion (Suara et al. 2015, 2016, 2017) by performing multi-device field studies using both Eulerian and Lagrangian measurements. The field studies showed that dynamics of Eprapah Creek is mainly controlled by the tide with a semi-diurnal pattern. Although observations are reliable, both Eulerian and Lagrangian devices cannot generally provide a complete picture of the environmental processes, especially those processes that vary spatially and temporally (Heemink et al. 2009). Therefore, the hydrodynamics of the estuary is required to be studied numerically with reasonable accuracy to ensure the efficient monitoring of environmental processes within the system.

Given the above discussion, this study focuses on evaluating the sensitivity of a 2D depth-averaged hydrodynamic modelling of a micro-tidal estuary to the bed roughness and bathymetry offset. Importantly, three scenarios using these parameters are used in the model calibration process to examine the impact of calibrating model using bathymetry offset. In these scenarios, the bathymetry adjustment and bed roughness were defined as uniform and constant for the domain of interest because of the small size of the domain. Numerical results are compared with time series of observed water level and velocity. As an important novelty in hydrodynamic modelling of shallow estuaries, high- and low-resolution Lagrangian drifter datasets are used to validate the model outputs in three experiments with different tidal conditions. This study has three main objectives:

• 1.

  Study the relative sensitivity of the model outputs to bed roughness and bathymetry offset. This analysis can assist not only to assess the effect of errors and approximation in bathymetry and bed roughness on model outputs, but also to provide an insight into the appropriate combination of parameter values.

• 2.

  Evaluating effectiveness of the simultaneous calibration of bed roughness and bathymetry offset against conventional calibration with a single parameter (e.g. bed roughness).

• 3.

  Exploring the capability of using Lagrangian drifter data, which is a new type of dataset in shallow estuaries, for the validation process.

The field of study was located at Eprapah Creek, which is a micro-tidal estuary in south east Queensland, Australia (longitudinal 153.293° E and latitude 27.574° S). This estuary consists of straight and meandering channels, and its morphology is characterized by irregular bathymetry with a maximum depth of 3–4 m (Figure 1). The estuarine zone is about 3.8 km long and due to the existence of mangroves, it is mostly sheltered from the direct wind. The mean annual rainfall of the region is around 1,179 mm and the cumulative rainfall for the last 7 days immediately prior to the experiment was approximately 30 mm at the weather stations within 20 km of the catchment (ABM 2020). Eprapah Creek discharges into the Moreton Bay in the Pacific Ocean, and the tide is the main driving force in this estuary.

A multi-device field study was conducted during the spring tide in July 2015. A full description of the field study and data analysis is presented in detail in Suara et al. (2017, 2018). The summary of instrumentation is given in Table 1.

Some key features of field data collection and other data used in this study are described below:

• Lagrangian measurements of velocity were obtained from deployment of the low- and high-resolution drifters in a relatively straight part of the estuary during the experiment (Figure 1). Low- and high-resolution drifters were sampled at 1 and 10 Hz with the position accuracy of 1–2 m and 2 cm, respectively (Suara et al. 2015). By using a finite difference approach, velocity data have been derived from the successive drifter positions. More details about the Lagrangian data collection and the method used to extract the drifter velocities are presented in the subsequent sub-sections. Velocity data from drifters were used to validate the hydrodynamic model.

• Eulerian measurements of water velocity were obtained from an acoustic Doppler current profiler (ADCP) and two acoustic Doppler velocimeters (ADVs) deployed upstream (approximately 2 km from mouth). A quality control was performed to remove the communication errors (i.e. data with correlation less than 60% and signal-to-noise ratio (SNR) smaller than 5 dB) in the ADV datasets. The ADCP datasets located in the air and depth sidelobe effects were removed. Moreover, the ADV and ADCP datasets were de-spiked by using the phase-spaced thresholding method. The spurious datasets identified by this method were less than 0.2% within the individual bin (Suara et al. 2018). The depth-averaged velocity extracted from the ADV by using the vertical-velocity curve is used to estimate the time series of discharge as the upstream boundary condition for the hydrodynamic modelling. Considering the uncertainties in the measured channel width, water level and velocity, the propagation of uncertainty to the time series of discharge is estimated to be around 2%. Depth-averaged water velocity obtained from the ADCP is later used for calibration of the model.

• The water level was collected manually using a staff, and its time series are used for the downstream boundary.

• Bathymetry was obtained from LiDAR data with a resolution of 5 m and a vertical accuracy of 0.3 m from Geoscience Australia. In addition, a survey of bathymetry was performed at different cross-sections of the estuary. The locations of all sites, where surveyed cross-section was carried out, as well as two graphs of measured cross-sections are shown in Figure 1.

During 48 h of field work, three experimental deployments of drifters were carried out. Two experiments (E1 and E2) were conducted in flood conditions with the tidal range of 1.75 and 2.25 m. The third experiment (E3) was performed during slack water with a tidal range of 1.7 m (each experiment period is shown in Figure 5). More information about the environmental conditions, including wind speed and direction and average water surface velocity, during the different experiments are presented in Table 2. Wind direction was measured clockwise from the positive streamwise direction downstream (shown in Figure 1(a)). To remove positioning errors at high frequencies where the SNR was less than 10, a low-pass filter with a cut-off frequency of 1.5 Hz was applied to the drifter data (Suara et al. 2015).

In this study, Delft3D Flexible Mesh (FM) is used to perform hydrodynamic modelling of the estuary. Delft3D-FM is an unstructured hydrodynamic model developed by Deltares. Using finite volume techniques on a staggered scheme, it solves the unsteady shallow water equations under a Boussinesq assumption in 2D depth-averaged or in three-dimensional modes, which are derived from three-dimensional Navier–Stokes equation (Deltares 2019). By using a flexible mesh, Delft3D-FM can be applied to complex geometries such as meandering estuaries. Delft3D-FM has been widely used in modelling flows in shallow water systems (Symonds et al. 2017; Thanh et al. 2017; Bomers et al. 2019).

In the present study, the hydrodynamic simulation is performed in a 2D depth-averaged mode. Performing a mesh independency experiment, a spatially variable unstructured mesh with an average 3.3 m resolution was found to be suitable for the model. This experiment was performed through evaluating the effect of the mesh size on the modelled velocity at several points. Further refinement of the mesh showed small changes in the value of the modelled velocity. Therefore, this mesh resolution was selected to have a balance between good representation of the system states and reasonable computation time. Higher-resolution grids were defined for the meandering parts of the estuary to capture the change in the direction of flow in these regions. Time-varying water level and discharge obtained from observations are used as downstream and upstream boundary conditions, respectively. The modelled domain is presented in Figure 1(b). The figure also shows the bathymetry relative to Australian Height Datum (AHD), which was interpolated onto the model grids via an interpolation method available in Delft3D-FM. Note that the contribution of the wind is assumed negligible since the tide is the main driving force during flood and ebb tides, although the high spatio-temporal variability of wind can become an important driver for the motion of the Lagrangian drifters in slack tides. A summary of the model set-up is presented in Table 3.

The main aim of this study is to analyse the sensitivity of a hydrodynamic model of a micro-tidal estuary to input parameters and eventually to minimize the systematic errors through calibration. A sensitivity analysis was performed by varying the bathymetry offset in the range between 0 and 1.2 m and Manning's coefficient from 0.016 to 0.050. These ranges of input parameters were selected based on physical characteristics of the system as well as the observed deviations between interpolated bathymetry at model grids and those manually obtained in the field work. In addition, the bathymetry offset and bed roughness with the values beyond the selected ranges in this case showed no improvement in the model performances. To provide a suitable amount of information for sensitivity analysis, 300 sets of two input parameters were generated with the step sizes of 0.002 for Manning's coefficient and 0.05–0.1 m for bathymetry offset. The step size for bathymetry has been selected since the sensitivity of the model to such changing in bathymetric data was obvious. In addition, the existence of curvature in the channel can cause an approximate increase of 0.002 in Manning's coefficient based on Chow (1959), and therefore, this value can be a reasonable step size for bed roughness for Eprapah Creek, which consists of straight and meandering parts. To have an insight into the approximate range of bed roughness, the Manning's coefficient can also be estimated through the volume balance equation (Amiri et al. 2016).

To evaluate the hydrodynamic model performance for representation of the essential physics of the system, the Pearson correlation coefficient (R) and root mean-squared error (RMSE) between model outputs and observations were calculated at each time step observations were available. These sensitivity indices were defined for the time series of water level (R_H and RMSE_H) and flow velocity (R_V and RMSE_V). In the case of RMSE_V, the magnitude of horizontal velocity is used. The flowchart of the sensitivity analysis is presented in Figure 2.

Using Lagrangian drifter data for the validation purpose is an important part of this study since this type of dataset is new in shallow estuaries and the number of studies which have used Lagrangian data to validate the hydrodynamic model of estuaries is very limited (Huhn et al. 2012; Mardani et al. 2020). Validation is performed by evaluating the modelled Eulerian flow velocity against the LR and HR Lagrangian drifter velocity data for three experiments with different tidal conditions (given in Table 2). The capability of the model to reproduce the flow in Eprapah Creek for different tidal conditions can therefore be evaluated.

It should be noted that the drifter velocity data are varying within a model grid cell since the time resolution of both low- and high-resolution drifters is higher than time resolution of the model, i.e. there is spatio-temporal difference between the model outputs and the drifter data. Therefore, the comparison between the model outputs and drifter velocity data is carried out through spatio-temporal averaging. In this application, there is a challenge in terms of what time resolution of observations compared with space and time resolutions of simulated data can meaningfully be used for the validation step. A simple analysis to address this issue is to determine how long it takes for a drifter to pass through a grid cell with its average velocity. Given this discussion, we use the time step of 200 s, which is 20 times larger than simulation time step, to compare the observed velocities with the simulated velocities.

The sensitivity of the model outputs to the input parameters is investigated based on the correlation coefficient and RMSE between time series of modelled and observed water level and flow velocity. Performing 300 simulation runs using different values of the bathymetry offsets and bed roughness described in the previous section, the contour plots of model performance indices are presented in Figure 4. These plots illustrate the range of the input parameters (dark red) in which the model produced the best results. Figure 4(a) and 4(c) show that the minimum values for RMSE_V and RMSE_H are 0.045 m/s and 0.08 m corresponding to Manning's coefficient of around 0.020 and bathymetry offset of 0.95 m. Water level is dependent on both bed roughness and bathymetry, while its sensitivity to the bathymetry offset is greater (Figure 4(c)). Similarly, both model input parameters have noticeable effect on the flow velocity, but the effect of bathymetry offset is more significant (Figure 4(a)). Hence, errors in bathymetric data impose larger uncertainties into the model outputs than do errors in bed roughness. This indicates that the adjustment of the bathymetry can significantly improve model performance indices. The optimum bathymetry offset based on both water level and flow velocity is 0.95 m. No further improvement in modelled water level could be achieved by using larger bathymetry offsets, while RMSE_V is still sensitive to the further variation of bathymetry offset. On the other hand, by comparing Figure 4(a) and 4(c), it can be noticed that errors in bed roughness introduce a larger uncertainty into the water level than flow velocity.

Similarly, R_V in Figure 4(b) converged to 0.98 at Δz_offset ≥0.9 m and small values of roughness. On the other hand, RMSE_V in Figure 4(a) reached a minimum value at Δz_offset = 0.95 m and small values of roughness. Further increase in bathymetry offset increased the deviation of the model velocity from the observation. This is because an increase in the bathymetry offset beyond the optimum level would result in an increase in the cross-sectional area and a reduction in the cross-sectional average velocity for a given discharge. In this case, therefore, the minimum value of the RMSE in the velocity is an important factor for selecting optimum model parameter in the calibration process. In other words, it is more efficient to perform the calibration of input parameters against flow velocity rather than water level since the sensitivity of flow velocity is greater than the sensitivity of water level to the input parameters.

Figures 4(b) and 4(d) show that the optimum ranges of the input parameters for R_V and R_H are much larger than those for RMSE_V and RMSE_H. Indeed, the linear correlation between model and observations can converge to its optimum for the certain values of bed roughness and bathymetry offset, yet RMSE needs more effort to converge to its optimum. Modelled and observed water level and velocity are highly correlated (R_V = 0.98 and R_H = 0.99) at a bathymetry adjustment of 0.9 m, and further changing in the bed roughness and bathymetry offset do not significantly affect the R_H and R_V. Both performance indices (RMSE and correlation coefficient) should be taken into consideration in sensitivity analysis in order to gain better insights into relative sensitivity of the model variables to the parameter variation, especially to determine which observed variable is more appropriate for the calibration of input parameters.

The errors associated with bathymetric data are strongly dependent on vertical accuracy and horizontal resolution of bathymetry survey. They have been quantified in a few studies, which are reported to be in the order of m (List et al. 1997; Van Der Wal & Pye 2003), and 1 m or more (Burningham & French 2011) depending on the technology and operational post-processing parameters. The optimum Δz_offset = 0.95 m obtained in this study is consistent with those in the literature. In addition, the comparison between manual bathymetry measurements and the LiDAR data for eight cross-sections in Eprapah Creek showed an offset ranging from 0.61 to 1.18 m. Therefore, the bathymetry offset obtained from the sensitivity analysis is representative of the dynamics of Eprapah Creek.

As the visualization of sensitivity analysis indicated, the sensitivity of flow velocity in terms of RMSE is larger than sensitivity of water level to the variation of input parameters. Therefore, the calibration of input parameters against the observed flow velocity is efficient in this study, and thus the cost function (GoF) was defined based on the modelled and observed velocity. The calibration was performed by using both the bathymetry offset and bed friction. However, to determine the influence of involving each input parameter in the calibration process, calibration with single parameters was tested as well. Table 4 presents the calibration statistics for each criterion. As it can be seen from this table, having only bed roughness as a calibration parameter (case C1) cannot significantly improve the model to represent the physical processes, and the deviation between modelled and observed velocities in terms of RMSE is noticeable (RMSE_V = 0.133 m/s). In this case, correlation between model outputs and observations is limited to R_H ≈ 0.90 and R_V ≈ 0.86. Compared with the calibration of bed roughness, calibration of bathymetry offset (case C2) resulted in noticeably smaller RMSE for both water level and flow velocity (RMSE_H = 0.117 m and RMSE_V = 0.048 m/s). The impact of bathymetry adjustment on RMSE_V is greater than its impact on RMSE_H, and after calibration, the model can satisfactorily represent the flow with a high correlation (R_V = 0.973) in the domain of interest. With inclusion of bed roughness to the calibration process (case C3), further improvement in model performance, especially for water level, can be obtained. By simultaneously calibrating bed roughness and bathymetry, the model can represent the water level and flow velocity with relatively small deviations from observations (RMSE_H = 0.078 m and RMSE_V = 0.045 m/s) and its output is highly correlated with observed data (R_H = 0.995 and R_V = 0.985).

Figure 5 shows the simulated and experimental time series and scatter plot of water level and flow velocity for the optimum input parameters (i.e. Δz_offset = 0.95 m and n = 0.021) obtained from calibration. The time series of water level show good agreement between model and observations in terms of amplitude and phase. Moreover, the scatter plot reveals that modelled and observed water levels are highly correlated (R_H = 0.99). Although there are some deviations between modelled and observed Eulerian velocity, Figure 5(b) and 5(d), and Table 4 indicate that the model reasonably captured the flow velocity. The deviation between the simulated and observed flow velocity is larger during flood tides, where the model underestimates the flow velocity, than ebb tide phase.

To have a better insight into model performance in different tidal phases, Table 5 presents the performance indices during ebb and flood tides. As shown in Table 5, the model has better performance in ebb tides for both water level and flow velocity with RMSE_H = 0.051 m and RMSE_H = 0.031 m/s, respectively. These differences could be related to the uncertainties in boundary conditions and the effect of wind forcing in field observations.

To validate the model, a comparison between Lagrangian drifter velocity and the Eulerian velocity of the model for three different experiments is performed. The scatter plot of modelled and drifter velocities for all experiments are illustrated in Figure 6(a) is related to HR drifter data and Figure 6(b) is for LR drifters. The deviation between Eulerian velocity and observed velocity in terms of RMSE is calculated. Generally, these results show an acceptable agreement between the model and both HR and LR drifter datasets. The drifter velocities are distributed on either side of the regression line for both LR and HR drifters which could be due to comparison of the depth-averaged velocities from model with drifter velocities which correspond to the velocity field near the water surface. Figure 6(a) shows that model mostly underestimated the velocity compared with HR drifter velocities, and it is generally expected in an open channel that the surface flow is larger than the depth-averaged velocity. However, Figure 6(b) reveals that generally the model overestimated the velocity compared with the LR drifter velocities. The other source of error that cannot be neglected arises from deriving the Lagrangian drifter velocities from their successive positions to represent a velocity value at a single point in the domain of interest. This error can be more significant for LR drifters than that for HR drifters due to relatively larger position errors of LR drifters (between 2 and 3 m) compared with the HR drifters (approximately 2 cm). The deviation between model and Lagrangian velocity obtained from HR drifters (RMSE_HR ≈ 0.070 m/s) is lower than that for LR drifters (RMSE_LR ≈ 0.082 m/s). Correspondingly, the Pearson correlation coefficient for HR drifters (R_HR = 0.87) is higher than that for LR drifters (R_LR ≈ 0.83). It is notable that the level of the correlation as well as the RMSE between model and drifter data varies with the phase of the tide in which the experiments were carried out. To clarify, the model performance metrics are calculated for each experiment (Table 6). It can be observed from Table 6 that during E1 and E2, the modelled velocity and drifter velocity are more highly correlated for the HR drifters. The correlation coefficient between model and LR drifters for E1 was high (with R_LR ≈ 0.87), while a lower correlation with the value of 0.71 was obtained for E2. These levels of correlation observed in the first two experiments are consistent with other research studies on comparison of drifter and modelled velocities. Huhn et al. (2012) reported that the correlation between hourly-averaged drifter and modelled velocities in a tidal driven estuary varies between 0.64 and 0.8, while Abascal et al. (2009) observed that this value range from 0.67 to 0.74 in their study on the coast of Spain.

Contrary to the first two flood experiments, the correlation coefficient for both LR and HR drifters reduces significantly during slack water E3 (with R_LR ≈ 0.30 and R_HR ≈ 0. 58). The model generally underestimated the flow velocity during this period. This significant reduction in regression can be explained by the combination of environmental factors during this experiment. As can be seen from Table 2, the average water surface velocity was small during slack, while wind speed was significant. Eprapah Creek is a semi-enclosed estuary, and it is surrounded by mangroves which caused a strong spatio-temporal variability of wind in this area. This variability of the wind can impose variation into the surface flow and, especially during slack, the effect of wind through wind-induced current on drifter motion could be noticeable (Suara et al. 2018).

Resonance, which is the reflection of the tidal signal between landmarks, has been reported to be significant in Eprapah Creek (Suara et al. 2019). Therefore, resonance can impose rapid fluctuations particularly in the direction of flow experienced by these drifters. Although the model might be able to capture resonance by imposing time series of discharge at the upstream boundary, there would be lower correlation between model outputs and Lagrangian drifter velocity at slack tide because wind as a driving force is not accounted for in this hydrodynamic model simulation.

The lack of correlation between modelled and observed data is not all due to the wind effect. It could also arise from position errors that are intrinsic to the GPS-tracked drifters, especially lower accuracy of the instrument at very low velocity like those experiments at slack tides. For example, the LR drifter data during slack water were dominated by noise since the magnitude of water surface velocity was less than 0.1 m/s, which is in the neighbourhood of the magnitude of inherent error in the speed estimate of the off-the-shelf LR GPS data-logger (Suara et al. 2018). Therefore, the inherent noise present in the drifter data also contributed to the differences between modelled and observed velocities. In this regard, limitations of drifters used in this study have been discussed in Suara et al. (2015).

In addition to the aforementioned factors, some other sources of uncertainties, such as imperfect knowledge about physics of the system, discretization errors, and initial and boundary condition errors, result in stochastic errors that can lead to divergence of the model outputs from observations. Among these uncertainties, modelling of the flow in the estuarine systems is dependent to a great extent on boundary conditions, and thus errors in boundary conditions can adversely affect model results (Wang et al. 2009). By employing DA, these uncertainties can be significantly reduced further. In summary, after calibration using both bathymetry and bed roughness, the model generally captured the flow velocity in the domain of interest. Importantly, this numerical study shows that the tidal force is the dominant driving force in Eprapah Creek, although effects of wind on water surface during slack water may be present.

In this study, Delft3D FM, as a hydrodynamic model, was employed to simulate the tidal flow of a micro-tidal estuary with complex geometry in a 2D depth-averaged mode. Inputs for hydrodynamic modelling were obtained from 48 h of field work in which both Eulerian and Lagrangian instruments were deployed and the model was forced by the tide. Bathymetry, used in the model, was extracted from LiDAR data; however, we observed noticeable deviations between them and eight manually surveyed cross-sections. A sensitivity analysis revealed that the errors in bathymetric data significantly affected the model outputs (i.e. water level and velocity). Performing a new bathymetry survey covering the entire domain was impractical due to cost and level of accuracy required; therefore, bathymetry errors and bed roughness were treated as calibration parameters to improve the model accuracy, while the manually surveyed cross-section was used to validate the bathymetry offset.

The sensitivity of the model outputs to the variation of bathymetry adjustment and bed roughness was performed. Visualization of the sensitivity analysis in terms of the Pearson correlation coefficient as well as RMSE for both water level and flow velocity indicated that the model performance is significantly sensitive to the bathymetry adjustment. Specifically, it was shown that errors in bathymetric data can result in model outputs that significantly deviate from observations. Furthermore, it was shown that the bed roughness and bathymetry can be efficiently calibrated against flow velocity rather than water level since flow velocity in this study is more sensitive, than water level, to variation of input parameters.

Results of the work showed that, considering bathymetry offset as the only calibration parameter led to significant improvement in model performance in terms of RMSE and the correlation coefficient between model outputs and observations. However, there were still differences between modelled and observed variables. Best performance, especially for water level, was achieved after simultaneous calibration of bathymetry and bed roughness. High correlations of around 0.98 and 0.99 were obtained for simulated and observed Eulerian velocity and water level, respectively. Furthermore, the results revealed that the model had a better performance during the ebb than the flood tides.

The model was validated using velocity data obtained from both low and high-resolution Lagrangian drifters. Validation was performed for three experiments with different environmental conditions. The results showed that performance of the model in comparison to the Lagrangian drifter velocity varied as a function of tidal phase. Generally, Eulerian velocity obtained from the model has good agreement with drifter velocity, especially during flood tides where tidal force was the dominant driving force. Moreover, model performance in comparison to high-resolution drifters was greater than that against low-resolution drifters. Lower correlation was obtained during slack tides, and this could be due to wind effect and low accuracy of drifters at low flow velocity.

Models calibrated through adjusting bathymetry can be used for real-life applications including operational and investigative purposes. Propagation of the tide can be improved by calibrating bathymetry, and thus such models can be used for storm surge forecasting (Zijl et al. 2013; Yan et al. 2015). For example, the hydrodynamic model of the Dutch Continental Shelf Model (DCSMv6) was calibrated by adjusting the bathymetry of the North Sea (Zijl et al. 2013). This model runs operationally to provide the water level forecasts every 6 h along the Dutch Coast (Zijl et al. 2013). Considering bathymetry errors as a calibration parameter can effectively improve the computed velocity fields (Cea & French 2012; Matte et al. 2017) and reduce the need to undertake a large-scale survey after events that can lead to a large bathymetry change in large domains.

This work provides, for the first time, a 2D hydrodynamic modelling of Eprapah Creek, which is calibrated and validated using both Eulerian and Lagrangian data. The outputs of this model, i.e. improved velocity fields, are currently being used for investigative purposes to drive the advection of materials. In particular, the velocity field is being used as an input to advection–diffusion models and for visualizing the structures responsible for Lagrangian transport in the system to gain better insights into material transport in this field. Moreover, this study is the first step towards the development of a comprehensive hydrodynamic model of this region which will be coupled with a DA system for a future operational forecasting system.

Model outputs are highly sensitive to bathymetry errors; therefore, conducting an accurate full-scale bathymetry survey is the ideal solution to mitigate this error. However, highly accurate bathymetry surveys are expensive and time-intensive, performing a new survey of bathymetry may not be practical for some environments. Unlike bathymetry, velocity measurements are usually easier to collect. A practical approach to estimating the bathymetry offset for an outdated or erroneous bathymetry is through observed flow velocity (Landon et al. 2014). Therefore, assimilation of the velocity measurements (depth-average and drifted-based data) can be employed to extract the correct bathymetry of the poorly gauged areas in semi-enclosed rivers and estuaries (e.g. Wilson & Özkan-Haller 2012; Landon et al. 2014; Almeida et al. 2018).

In this study, the bathymetry adjustment and bed roughness were defined uniform and constant for the domain of interest because of the small size of the domain. Further improvement in the calibration process, particularly for larger systems, could be achieved by considering spatially varying input parameters. Second, although the 2D depth-averaged hydrodynamic model reasonably represented the physical processes in Eprapah Creek, performing 3D modelling can bring more insights into vertical mixing dynamics of this estuary. The interaction of different environmental forcings, such as wind, with tide was not studied in this work, and thus better representation of the flow can be achieved by considering wind in the hydrodynamic modelling especially during slack water. Having an established model with low systematic error via the calibration process, the main objective of the future study is to effectively mitigate random errors by performing DA using both Lagrangian and Eulerian observations as well as further modelling for flow forecast purposes.

Discussion and advice from Prof. Hubert Chanson and Mr Firmijn Zijl are gratefully acknowledged. The authors would also like to thank the editor and three anonymous reviewers for their constructive comments on the manuscript. This project is supported through Australia Research Council Linkage grant LP1501072 and Discovery Project grant DP190103379.

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