COMPARISON BETWEEN MIKE 21 FM, DELFT3D AND DELFT3D FM FLOW MODELS OF WESTERN PORT BAY, AUSTRALIA

The performance of three different hydrodynamic modelling packages is compared in this study, namely Delft3D, Delft3D FM (both developed by Deltares) and MIKE 21 FM (developed by DHI). Delft3D and MIKE 21 FM are internationally known software packages while Delft3D FM (formerly known as D-Flow FM) is a relatively new package. The models use structured approaches (Delft3D), unstructured approaches utilising triangular and quadrilateral elements (MIKE 21 FM) and unstructured approaches utilising elements ranging from linear to six sided (Delft3D FM). Models of Western Port, Australia, were developed using the three different packages to allow a comparison of performance and to determine if there are any differences in using structured versus unstructured approaches. Model performance has been assessed based on model calibration, representation of channel flows and computational efficiency. Despite the inherent differences in the grid configuration and the implementation of the numerical schemes between structured and unstructured approaches, both approaches have been shown to be able to accurately predict hydrodynamic conditions in a complex estuarine environment. The unstructured approach was found to be the most computationally efficient both when run on multiple cores (MIKE 21 FM was the most efficient) and when run on a single core (Delft3D FM was the most efficient).


INTRODUCTION
There is a wide range of numerical modelling software available for marine, coastal and estuarine modelling.Many of the models have been developed by universities and research institutes, and models have also been developed by private companies.Some of these models have subsequently been made publically available and used in commercial projects, while some models are used primarily for research within academia.The software considered in this study is all publically available and can be adopted as part of commercial projects.
The models tend to be based on similar mathematical concepts.For example, most hydrodynamic models utilise the Navier-Stokes equations.The models can either adopt a structured or unstructured grid, with both approaches being widely adopted internationally for a range of applications.The structured grids can be rectilinear or curvilinear and typically adopt a finite difference solution scheme, while the unstructured grids have a flexible mesh (mesh elements can have a variable shape and size over the model domain) and typically adopt a finite volume solution scheme.
The difference between structured and unstructured approaches is not well documented in the existing literature.The relative difference is dependent on the particular grid/mesh configuration and has to be assessed on a case by case basis.This paper considers a case study where models have been developed of Western Port, a large estuary in Australia using three different packages.The models have been configured to have comparable resolutions, similar offshore extents and they are driven by the same offshore boundary conditions, to allow comparison of performance and determine if there are any differences in using structured (Delft3D) versus unstructured (Delft3D FM and MIKE 21 FM) approaches.Model performance has been assessed based on the model calibration, representation of channel flows and computational efficiency.

GEOGRAPHICAL SETTING
Western Port is located approximately 60 km south-east of Melbourne in Victoria, Australia.It is an estuary of irregular shape extending about 40 km east-west and 40 km north-south overall and covering an area of about 680 km2 (with about 270 km 2 of intertidal area exposed at low tide), with two ocean entrances.Phillip Island separates the two ocean entrances, and French Island divides the estuary into a North Arm and East Arm.
The geographical setting of Western Port is provided in Figure 1, along with discretisation into segments based on physical characteristics as derived from Marsden et al (1979).The boundaries of the segments are generally situated at relatively narrower sections of the waterway, with four such constrictions identified in Figure 1.
The configuration of Western Port with two entrances, two large islands within the bay and extensive areas of shallow intertidal flats results in the bay being subject to complex physical processes.In the eastern part of the Upper North Arm the tidal channels converge and are separated at low water by a 'tidal divide' (see red dashed line in Figure 1), which plays an important role in the tidal flows elsewhere in the bay (Marsden and Mallett 1975;Miles 1976;Rosengren 1984).The complexity of the processes and the configuration of the bay makes it a challenging area to accurately represent in a numerical model and therefore a suitable location to compare different numerical modelling approaches.

MODEL DESCRIPTION
The three hydrodynamic models were setup as two-dimensional (2D), depth-averaged models.In a depth-averaged 2D configuration, all three models are based on the numerical solution of the depth-integrated incompressible Reynolds averaged Navier-Stokes equations.The models consist of continuity, momentum, temperature, salinity and density equations.A Boussinesq assumption is applied in all the models, in which it is assumed that momentum transfer caused by turbulent eddies can be modelled with an eddy viscosity.
The spatial discretization of the equations in the MIKE 21 FM and Delft3D FM models, which adopt unstructured grids, is performed using a cell-centred finite volume method (DHI 2014a; Deltares 2015).For these models the time integration of the shallow water equations is performed using an explicit scheme (Delft3D FM adopts a combined implicit and explicit scheme, for further details refer to Kernkamp et al. (2011)).Due to the stability restriction using an explicit scheme, a dynamic time step is adopted in both models whereby the model specifies the time step for each computational step to help maintain model stability and performance.In the horizontal plane MIKE 21 FM uses an unstructured mesh comprising triangles or quadrilateral elements, while Delft3D FM uses an unstructured mesh comprising linear (one-dimensional), triangular, quadrilateral, and polygonal cells with at most 6 sides.
The Delft3D model, which adopts a structured grid, performs the spatial discretization of the equations using a cell-centred finite difference method (Deltares 2014).In this model the time integration of the shallow water equation is performed using an implicit scheme.Implicit schemes are unconditionally stable (although non-linearities in the governing equations can limit stability), allowing a larger Courant number than explicit schemes which in turn allows a larger time step.As such, a fixed time step is applied in the Delft3D model which is specified by the modeler and is typically defined based on the courant number.In the horizontal plane a structured grid is used comprising rectilinear or curvilinear cells.

Bathymetric Data
It is preferable to use the most recent and highest resolution bathymetric maps in a numerical model to ensure the model can provide as realistic a representation of the actual environment as possible.For this comparative assessment of different numerical models it is important to adopt the same bathymetric data in each model.
High resolution bathymetric and topographic data was available throughout most of Western Port from airborne LiDAR and multi-beam echo sounder surveys.The extent of this combined data set is provided in Figure 2. Limited recently measured bathymetric data was available for the area offshore of Western Port and so the MIKE C-MAP data was adopted for this area.The different bathymetric data sets utilised were corrected to Australian Height Datum (AHD) and then used to create a Digital Elevation Model (DEM) which could be interpolated onto the model grids/meshes of the three models.

Offshore Boundary
Sensitivity testing of various offshore boundary locations and configurations was undertaken to optimise the offshore boundary for the models.The testing found that the optimum boundary configuration for the models was located approximately 30km to the west of the Western Entrance to Western Port, 30km to the south-east of the Eastern Entrance to Western Port, and 35km offshore of Phillip Island (Figure 3).As such, all three models were set up with this same offshore boundary location.
The only long term water level data available within Western Port or in the immediate offshore surrounds has been collected at Stony Point (Figure 1).The lack of measured water level data close to the offshore model boundaries meant that the best approach to define water level boundary conditions was to adopt predicted water levels from a Global Tidal Model.The astronomical tide water level boundary condition adopted was based on the 0.125° resolution Global Tidal Model developed by DTU Space (DHI 2014b).This model contains 10 tidal constituents including the M4 shallow water constituent and has been used to represent the astronomical tidal boundaries for all three of the models.

Model Domains
To accurately represent the important hydrodynamic processes in Western Port it was determined that the resolution of the model grid/mesh domains should be as follows:  Western Entrance Segment, Confluence Zone and Lower North Arm need to be high resolution to ensure that the channel cross-sectional area is accurately represented. Upper North Arm and East Arm also need to be relatively high resolution to be able to represent the wetting and drying of the intertidal area and to ensure the channel cross-sectional areas are accurately represented.Sufficient resolution is also required to realistically model the flow through the relatively narrow (400 m) eastern entrance to the ocean. The area offshore of Western Port could be relatively coarse, although near the entrance channels the resolution must increase to accurately replicate both the complex bathymetry and hydrodynamic processes in these areas.The three models were setup based on the required model resolutions detailed above and with the aim that they should have the same extent and comparable grid resolutions wherever possible.A summary of the range of grid cell or element arc lengths in the models are provided in Table 1.Further details of the computational grid/mesh configurations for each model are provided in the following sections and in Figure 4. Delft3D.Given the configuration of Western Port with wide/deep channels, two entrances and large intertidal areas, a curvilinear grid technique in combination with domain decomposition (DD) was adopted.DD is similar to grid 'nesting', but with improved computational efficiency as there is no overlapping of the different resolution grids.The curvilinear-DD grid technique allowed the total number of grid cells to be minimised while:  Defining high resolution in the areas of interest (entrance channels and main channels). Defining low resolution at the outer limits or in less relevant areas of the model. Providing dynamic coupling between model areas (domains). Providing the opportunity to further refine areas of interest as required. Typically allowing the grid orientation to follow the main channels.The three DD boundaries between grid 1 and grid 2 were required due to grid development restrictions in Delft3D rather than to allow a change in model grid resolution.
MIKE 21 FM.The MIKE 21 FM mesh was constructed with triangular elements, as this is the preferred flexible mesh construction approach for this software.The mesh was set up with coarser elements offshore of Western Port and then increased resolution within Western Port and specifically higher resolution in the main channels.
Figure 4 shows how the resolution of the mesh varies, with finer elements in the main channel and coarser elements on the adjacent shallower subtidal and intertidal areas.When compared to the structured Delft3D model grid it can be seen how the adoption of an unstructured grid allows increased variability in the spatial grid resolution relative to a structured grid.
Delft3D FM.Delft3D FM is capable of combining the use of curvilinear cells in deeper channels and triangles in topographically complex areas (such as intertidal areas).In deeper channels, high flow velocities can occur with uniform flow directions, which can be modelled more efficiently with curvilinear cells compared to a triangular grid system.Intertidal areas (which comprise channels, creeks and intertidal flats) can be modelled using the curvilinear grid technique for the channels and triangular cells for the flats.This technique has the benefit of improving model accuracy and efficiency in these areas.
The Delft3D FM mesh was designed with a similar curvilinear system for the main channels as in the Delft3D grid and with triangular elements in areas adjacent to the main channels (Figure 4).Delft3D FM has the advantage over Delft3D in that it can have variable resolutions in one model domain.This is beneficial to prevent overly high resolution in less relevant areas and to therefore reduce computational time.
Summary.The MIKE 21 FM model had approximately 100,000 computational elements, Delft3D had 88,000 and Delft3D FM had 80,000.The MIKE 21 FM model had more elements than the Delft3D and Delft3D FM models, because adopting a curvilinear grid with rectangular elements allowed variable resolution in the X and Y directions, but this was not possible in MIKE 21 FM where triangular elements were used.The variable resolution in the X and Y directions can enable a higher grid resolution to be adopted across a channel (where the variability in bathymetry and flow conditions is typically greater) and a lower grid resolution in the longitudinal direction of the channel (where the variability in bathymetry and flow conditions is typically smaller).This has the potential of improving the efficiency of a model without reducing the accuracy.

Bed Roughness
The spatially varying bed roughness coefficient is typically adjusted as part of the model calibration process, with consideration of the bed properties of different areas with varying bed sediments, bedforms (ripples, sand waves, etc) and vegetation (seagrass, mangroves, etc).For the purpose of this comparative assessment a spatially constant bed roughness was adopted.A range of constant bed roughness values were tested using an iterative approach and the value which gave the best calibration was adopted.Based on this, a Manning's n roughness coefficient of 0.025 was adopted for all three of the models.

Calibration
It is important to determine whether the models are capable of accurately representing the hydrodynamic conditions which occur in Western Port.Therefore, the three 2D hydrodynamic models were calibrated using:  current data collected using two Acoustic Doppler Current Profiler (ADCP) devices located in the main channel of the Lower North Arm; and  measured water level data from Stony Point.The locations of the two ADCP sites and the location of the Stony Point tide gauge are provided in Figure 5.
A harmonic analysis, considering up to 69 constituents and including all shallow water constituents, was undertaken on the measured data to remove any residual water levels and currents to provide results resulting from the forcing of just the astronomical tide.This allowed a direct comparison with the modelled results as the models were run with only an astronomical tidal forcing, which is considered to be acceptable for model comparison purposes.
Water Levels.Modelled and measured water levels at Stony Point are compared for the three different models in Figure 6.In addition, calibration statistics comparing the model results to the measured water levels for the three models are provided in Table 2.
All three of the models show comparable levels of calibration relative to the measured water level data, with a generally good agreement over the 29 day simulation.The differences in high and low water levels were less than 4% of the tidal range for all of the models, with the smallest differences for the Delft3D model (less than 2%).All of the models tended to slightly over-predict the high water levels while slightly under-predicting the low water levels (that is, being lower), meaning that they also tended to slightly over-predict the tidal range.The phasing of the tide was good for all the models, with a mean high water phase difference over the 29 day simulation of less than 5 minutes for both high and low waters.
There were periods when all of the models consistently over or under-predicted the water levels, namely:  the high water levels were slightly over-predicted by the models during the transition from spring tides to neap tides;  for one of the spring-neap cycles in the lunar cycle the high water levels were under-predicted when transitioning from neap tides to spring tides; and  for one of the neap cycles in the lunar cycle the high water levels were over-predicted.As these differences were consistent in all of the models this suggests that they are errors in the boundary conditions (astronomical components) rather than errors associated with the numerical models.Tidal Currents.Modelled and measured tidal currents at the northern ADCP site are compared for the three different modelling systems in Figure 7.In addition, calibration statistics for the three different modelling systems are provided in Table 3.The plots and calibration statistics demonstrate that all three of the models have comparable levels of calibration relative to the measured tidal current speed and direction data, with a generally good agreement over the 29 day simulations.
The mean peak flood and ebb differences were less than 0.02 m/s for all of the models, with the smallest differences using Delft3D.The root-mean-square (RMS) errors, along with Figure 7, show that for all three models there was some variability in the current speed calibration through the 29 day simulations.The periods with larger differences between the modelled and measured current speeds occurred during neap tides when the model over-predicted current speeds.These differences coincided with the periods when the model tended to over-predict high water levels, indicating that this difference is again at least partially related to the model boundary conditions.All of the models showed a good calibration for tidal current directions, with differences of less than 3° for both peak flood and ebb directions.A similar level of calibration was achieved for all three of the models at the southern ADCP site.

Hydrodynamics
Comparison of the peak tidal currents in the Lower North Arm predicted by the models shows that the speeds and spatial pattern predicted by the models are similar (Figure 8 to Figure 10).There are some differences between the flows in the MIKE 21 FM model relative to the flows in the Delft3D and Delft3D FM models, for example the MIKE 21 FM model predicts lower current speeds and a wider flow in the channel to the east of the main channel compared to the other two models.The differences between the model grid types and configurations will result in some differences in how the models represent the bathymetry within Western Port.The differences in the flow noted above can be attributed to the different representation of the bathymetry between the triangular elements used in MIKE 21 FM and the curvilinear grid adopted in this area in the other two models.
It is important for numerical models to be able to accurately represent how intertidal areas flood and drain through a tidal cycle.The wetting and drying of the intertidal areas can influence both the local hydrodynamics and can be an important process for both the deposition and resuspension of fine grained sediment on the intertidal flats.A comparison of wetting and drying in the three models was undertaken and it was found that the models showed similar intertidal depth and current patterns over a tidal cycle.At low water the extent of the intertidal zone which was exposed was comparable between the three models.It was therefore concluded that despite differences in the schemes the models adopt to represent flooding and drying, the models all predict similar flooding and drying of the intertidal areas in Western Port.

Computational Efficiency
For complex models which can take many hours to run it is important to consider their computational efficiency as if the simulations are impractically long the modelling approach may not be realistic.This is especially the case if additional more numerically demanding modules such as sediment transport and morphology will be required as part of the study.
A structured approach has a more efficient model solution than an unstructured approach but this can be counteracted by the unstructured approach having an explicit numerical scheme with a dynamic time step5 , while the structured approach has an implicit numerical scheme with a fixed time step.In addition, an unstructured approach provides more flexibility when constructing a complex model domain compared to a structured approach using square or rectangular cells.As such, the relative difference in run times between a structured or unstructured approach would depend on the particular grid/mesh configuration and therefore has to be assessed on a case by case basis.
All of the model domains developed as part of this study were optimised in terms of model efficiency to reduce the model run times as much as possible.The three different models were run for a 31 day period on a single core computer and multiple core computers6 ; with results from tests shown in Table 4 7 .
Models which adopt unstructured approaches including MIKE 21 FM and Delft3D FM can make use of machines with multiple cores by running in either Open Multi-Processing (OpenMP) or Message Passing Interface (MPI) modes.OpenMP is where the processing for the computation is split between the different cores, and has been used for the multiple core model simulations outlined in Table 4.This approach generally gives a good scalability for 2-4 cores but diminishing scalability for a larger number of cores (Sørensen et al. 2010).MPI is where the domain is partitioned into multiple equal size areas and each partition is run on a separate core.The scalability of this approach is dependent on the complexity and uniformity of the model mesh, but it can give close to a linear speedup factor as the number of cores is increased.This approach was tested for the Western Port MIKE 21 FM model, and it was found that the MPI model ran in approximately half the time of the OpenMP model when 16 cores were utilised8 .The Delft3D model had to be run in MPI mode as it cannot utilise the OpenMP approach, this allowed each of the DD domains to utilize a single core meaning that the model ran using two cores.
Note that the model run time results in Table 4 are not directly comparable between models as the models have varying numbers of grid cells/elements and so the model efficiency is also presented whereby the run time has been divided by the total number of grid cells/elements.The MIKE 21 FM model had the most computational elements with 100,000, Delft3D had 88,000 and Delft3D FM had 80,000.The model efficiency shows that Delft3D FM was the most efficient model when run on a single core and the MIKE 21 FM model was the most efficient model when run on multiple cores.The results highlight the significant improvement in model efficiency when the MIKE 21 FM model is run on multiple cores compared to on a single core.It is also important to note that in terms of computational efficiency based on the number of cores utilised, Delft3D performed substantially better than the other two models as the 12 hours computational time was with only two cores.As such, on a multiple core computer, Delft3D would be able to have multiple simultaneous runs (for example, on an eight core computer, four simulations could be completed in about the same time that one simulation could be completed using MIKE 21 FM or Delft3D FM).

FINDINGS
This case study has found that all three of the numerical models considered predict similar hydrodynamic conditions within Western Port which also agrees well with measured data.As such, despite the inherent differences in the grid configuration (and therefore the bathymetry) and the implementation of the numerical schemes between structured and unstructured models, both approaches have been shown to be able to accurately predict hydrodynamic conditions in a complex estuarine environment.Therefore, with due consideration of the bathymetry and key processes, a modeller should be able to develop an adequate grid/mesh using either structured or unstructured approaches, assuming that the model used is suitable for simulating the processes of interest.
The unstructured models were found to be the most computationally efficient models both when run on multiple cores (MIKE 21 FM was the most efficient) and when run on a single core (Delft3D FM was the most efficient).However, for the multiple core testing the structured model could only utilise two cores while the unstructured models utilised eight cores.Therefore, on a multiple core computer (e.g.eight cores) the structured model (Delft3D) could be able to have multiple simultaneous runs (four runs each using two cores) which would complete in a similar timeframe to a single simulation with the unstructured models utilising all eight cores.

Figure 1 .
Figure 1.Geographical setting of Western Port, Australia.

Figure 2 .
Figure 2. Extent of the high resolution aerial LiDAR and multi-beam survey around Western Port.

Figure 3 .
Figure 3. Extent and configuration of the Delft3D hydrodynamic model grids.

Figure 4 .
Figure 4. Close up of the Delft3D (left), MIKE21 (Middle) and Delft3D FM (right) model grids in the Lower North Arm.

Figure 5 .
Figure 5. Location of measured data used for the model calibration.

Figure 6 .
Figure 6.Modelled (all three models) and measured water levels at Stony Point for 29 day simulation, 7 day neap period (middle) and 7 day spring period (bottom).

Figure 7 .
Figure 7. Modelled (MIKE 21, Delft3D and D-Flow) and measured tidal currents at northern ADCP for 7 day spring period (current speed at top, and current direction at bottom).

Figure 8 .
Figure 8. MIKE 21 FM modelled tidal currents at peak ebb during a large spring tide in the Lower North Arm.

Figure 9 .
Figure 9. Delft3D modelled tidal currents at peak ebb during a large spring tide in the Lower North Arm (the white line is the DD boundary location, where due to the interpolation scheme adopted a small gap appears).

Figure 10 .
Figure 10.Delft3D FM modelled tidal currents at peak ebb during a large spring tide in the Lower North Arm.