1.    Background

A high-resolution 3D operational hydrodynamic ocean model around the entire Australian coastline, or national re-analysis, does not currently exist. To date, coastal case studies have either been performed on BLUElink outputs directly, or downscaled models have been developed which are nested within BLUElink. These latter models often require numerous nests to achieve acceptable nesting ratios when downscaling from the global model (10 km resolution) to the area of interest at the coast (100 s metres). A high resolution coastal model (of order 1 km resolution) would allow these ‘bridging’ nests to be removed when downscaling, and would also provide a suitably coastally optimized platform (which BLUElink is not) within which to conduct coastal oceanographic assessments.

Although an operational national modelling platform does not exist, there have been several attempts to develop such models. The technical barriers to a national model are many, and many of the current attempts do not adequately address them. The size of the domain is the dominant obstacle, i.e. how to achieve high coastal resolution without an impossible computational burden. A model that encompasses the rectangular size of the Australian continent is at a minimum approximately 5000 x 4000 km. If this were to be resolved at 1 km resolution a grid of 20 million surface cells results, which is a massive computational burden (most coastal applications being less than 100,000). Additionally, large computational inefficiencies result due to the large amount of land in a rectangular domain that must be integrated over, and there exists areas which are well beyond the coastal margins of Australia and not of interest, but must be included in a rectangular domain. It is advantageous to align cells with the principal directions of flow (alongshore and cross-shore), and this is not possible everywhere in the domain using rectangular meshes.

An alternative is to use a cyclic polar grid that wraps around the Australian continental shelves like a ‘ribbon’. This approach naturally excludes any land cells, and all cells in the domain are within the coastal margin and therefore of interest. This approach uses significantly less cells, where all cells from the rectangular approach that are over land or the in areas not of interest can be pushed into the coastal margins. This approach potentially could achieve a target resolution of 1 km with significantly less compute resources than a rectangular grid, however, the generation of these grids is extremely complex, time consuming, and ultimately does not possess the control to place required resolution in desired regions.

An unstructured model does not suffer these technical limitations. The unstructured mesh can be fitted to the coastline, with any arbitrary limit offshore, hence excluding all land cells. Resolution, and resolution transition, can be exactly placed where desired via the use of a weighting function. This function may itself be a function of bathymetry, distance from the coast, tidal amplitude, population density, eddy kinetic energy, any combination of these, or indeed any arbitrary metric. The unstructured approach is particularly suited to a national modelling exercise, and could practicably produce a mesh with coastal resolution of ~1 km. This project aims to produce such a model.

The approach of the National Modelling Initiative aims to produce similar products to the eReefs project, specifically a near real-time operational model (capable of short term forecasts if adequate forcing data is available) which appends daily output to a growing historical archive, and supported by a data assimilating re-analysis of that archive. The national model is proposed to be subdivided into a number of discrete ‘tiles’ representing distinct geographical regions (e.g. Great Barrier Reef, East Australian Current, Tasmania etc.) which are then assembled into a national model using 2-way nesting. Each tile could be operated in isolation for calibration and validation purposes, thus increasing throughput to result in a more constrained calibration, and a calibration that is region specific.

2.    v2 pilot model

The national model developed within this project uses the unstructured model COMPAS (Coastal Ocean Marine Prediction Across Scales; Herzfeld et. al., 2020,

https://research.csiro.au/cem/software/ems/hydro/unstructured-compas/ ). COMPAS is a three-dimensional coastal ocean model designed to be used at scales ranging from estuaries to regional ocean domains. The equations of motion are discretized on arbitrary polygonal meshes according to the dynamical core used by Ringler (2013) in the global MPAS model (Model Prediction Across Scales). This framework is described in Thuburn et al. (2009) and Ringler et al. (2010) and is hereafter referred to as TRiSK (Thuburn Ringler Skamarock Klemp). The TRiSK formulation is a generalization of the standard Arakawa C-grid scheme to unstructured grids.

COMPAS is based on the framework outlined in Herzfeld (2006); this model is a finite difference model that uses an underlying unstructured coordinate system. The generalization of this unstructured system greatly simplifies the transition from finite differences to TRiSK numerics. As such, the horizontal terms in the governing equations (momentum advection, horizontal mixing and Coriolis) are discretised using the TRiSK numerics, whereas the pressure gradient and vertical mixing are discretised using the finite difference approach outlined in Herzfeld (2006).  COMPAS therefore inherits much of the functionality of the host model of Herzfeld (2006), including numerous turbulence closure schemes (k-e, k-w, Mellor-Yamada 2.0 and 2.5) and open boundary infrastructure. COMPAS is available open source at https://github.com/csiro-coasts/EMS, with documentation provided at https://research.csiro.au/cem/software/ems/ems-documentation/ .

Meshes used with COMPAS are required to be a set of polygonal cells which are orthogonal (cell edges must be perpendicular to the line joining adjacent cell centres), centroidal (triangle vertices must be positioned at the centre-of-mass of each polygonal cell) and well-centred (cell vertices must lie in the interior of their associated dual triangles). These polygons are typically constructed from Voronoi diagrams, of which the dual is a Delaunay triangulation. Additionally, meshes must also conform to the geometry of coastal and open boundaries, and adhere to user-defined rules (a weighting function) that specify the mesh resolution. The generation of grids conforming to these requirements is not trivial, requiring the construction of highly optimised Centroidal Voronoi Tessellations (CVT’s) and ‘non-obtuse’ Delaunay triangulations. COMPAS uses the unstructured meshing library JIGSAW (Engwirda, 2017; https://github.com/dengwirda/jigsaw) to accomplish these tasks, resulting in high quality meshes that support the requirements of the TRiSK numerics.  JIGSAW is incorporated as an ‘inline’ component in the COMPAS framework, facilitating automated domain generation and model configuration. Various weighting functions are available inline in COMPAS; in this project we use a distance from the coast weighting. Initially the target coastal resolution is 2 km, with the view that this will be increased to 1 km when any emergent model development issues are streamlined. This 2 km application is referred to as the version 2 (v2) model.

The coastal weighting is considered the most appropriate if the objective of the national model is to provide system characterization of coastal processes, and provide downscaling forcing for coastal models at the estuarine or embayment scale. It may be necessary to refine resolution over areas not near the coast (e.g. reefs in the GBR). This is possible in COMPAS, as inline weighting also includes distance from a user defined polygon, or hybrid polygonal and coastal distance weighting.

The resolution map of the 2 km model is shown in Fig. 2.1. While this mesh maintains resolution of ~2 km at the coast, it decreases to as much as 25 km at the open boundaries. The islands in the Timor and Arafura Sea regions naturally produce high resolution at these northern open boundaries (due to coastal proximity), and this must be over-ridden so as coarser resolution is generated to maintain acceptable nesting ratios for boundary forcing.

Figure 2.1. Resolution of the National Model mesh.

The model uses bathymetry from the Geosciences Australia (2002) database, with regions outside this extent filled using the dbdb2 (Naval Research Laboratory Digital Bathymetry Data Base https://www7320.nrlssc.navy.mil/DBDB2_WWW/) global database. These are supplemented with the Beaman (2010) 100 m dataset within the Great Barrier Reef, and a 30 m northern Australia dataset (https://ecat.ga.gov.au/geonetwork/srv/eng/catalog.search#/metadata/121620). The GBR and NW Australia bathymetry datasets were interpolated onto the mesh using a bi-linear interpolation, and the remainder of the mesh was interpolated using the combined GA and dbdb2 dataset using a natural neighbours non-Sibson method. This allowed for accurate bathymetric values in areas requiring extrapolation (e.g. where 4 data points did not surround a model mesh point as required by bi-linear interpolation).

The model is nested within the global OFAM model (Oke et al., 2008) which supplies currents, sea level, temperature and salinity on the open boundaries. Low frequency sea-level from the global model is superimposed with 8 tidal constituents (M2 S2 N2 K2 K1 O1 P1 Q1) from the TPXO9v1 1/6° global model (Egbert and Erofeeva, 2002; http://volkov.oce.orst.edu/tides/otps.html) and applied at the boundary of the regional model. The boundary condition used is that of Herzfeld and Andrewartha (2012), modified for unstructured meshes. Surface fluxes are derived from the Australian Bureau of Meteorology’s operational atmospheric models (ACCESS-A; Puri et. al., 2012, http://www.bom.gov.au/nwp/doc/access/NWPData.shtml) at a resolution of 12 km. Heat fluxes applied at the surface boundary condition were computed from standard meteorological variables provided by ACCESS-A (wet and dry bulb temperature, air pressure, wind speed and cloud amount) using short and longwave calculations outlined in Zillman (1972) and the bulk method for sensible and latent heat using bulk coefficients of Kondo (1975). For the surface freshwater fluxes, precipitation was provided by ACCESS-A and evaporation was computed from the latent heat flux. Wind speed was converted to stress using the bulk scheme of Large and Pond (1981).

The model uses 53 layers in the vertical with 1 m surface resolution, the k-e turbulence closure model of Burchard et al, (1998) and an unstructured implementation of the ULTIMATE QUICKEST advection scheme for tracers. A bi-harmonic scheme was used for horizontal viscosity. Regionalized viscosity and diffusivity were applied, with generally more mixing supplied in the energetic regions of north west and east Australia.

The model is currently running in near real-time, with an archive created dating to September 2020. Output is generated in the Climate Forecast (CF) compliant netCDF standard for unstructured models, UGRID 1.0 (http://ugrid-conventions.github.io/ugrid-conventions/ ), which may be visualized using GODIVA 3 if data is hosted on OpenDAP.

3.    Outputs

The model is currently operating in near real-time beginning Sep 2020 on HPSC infrastructure, using 72 processors which delivers a runtime ratio of 2.85:1. Surface temperature with 3D currents and surface salinity with 2D currents are shown in Fig. 3.1 for Mar 2017. The distributions show no spurious behaviour, e.g. checkerboard patterns indicating instability or unrealistically large currents.

Figure 3.1. (a) Surface temperature and currents and (b) surface salinity and depth averaged flow for 12 Mar. 2017.

The model is compared with the eReefs GBR4 model for 5 Mar 2017 in Fig. 3.2 and 3.3 (note the different time zones). The broadscale patterns of flow in Fig. 3.2 are comparable, however, the GBR4 model exhibits a stronger North Queensland Current and more structure offshore in the Coral Sea. This is primarily due to the higher resolution (~4 km) over the whole domain in GBR4, whereas the National Model only has higher resolution within the coastal band and is resolved at ~25 km over much of the Coral Sea. The parameter configurations differ between the two models also; notably GBR4 is relaxed to global model density field in deep water to preserve currents driven by dynamic height gradients, whereas this relaxation is currently absent in the National Model. It is these parameter optimisations that require much attention during a calibration phase of the National Model.

The GBR4 flow is more chaotic at the NSW border, with the presence of an eddy in the SE corner, whereas the EAC off Fraser Island is stronger in the National Model.  Some boundary specification error likely contributes to GBR4 in this region, since the EAC exits the domain through the southern cross-shelf boundary and is also associated with a mesoscale eddy field. It is difficult to perfectly specify boundary forcing data in these energetic regions, and some over-specification error can result. Being of national scope, the National Model does not need to prescribe these cross shelf open boundaries, and this area demonstrates a clear benefit of modelling the entire shelves continuously. GBR4 must also use an open boundary across Torres St (which again is absent in the National Model), and to a lesser degree we see a deviation in flow patterns in this area between the two models.

Figure 3.2. Surface currents for the national model (top) and eReefs GBR4 (bottom) for 12 Mar. 2017. Note the different times due to time-zones used. Left column is eastward current, middle is northward current and right are the flow vectors.

Fig. 3.3 compares sea level, salinity and temperature between the two models. Sea level is again comparable, with an ebb tide seen in Broad Sound in both models. The peak tide is weaker in this region in the National Model, also the background sea levels are offset at the offshore open boundary. This indicates some tuning of the tidal forcing is required in the National Model. Differences at the southern and Torres St open boundaries are also observed, with lower sea level associated with the eddy near the SE boundary in GBR4 and weaker ebb tide through Torres St.

Salinity is broadly consistent between the two models, however, the coast is fresher in some regions in GBR4 due to freshwater river inputs which are yet absent in the National Model. This is most notable at the mouth of the Fly River, where low salinity in the National Model is retained from the initial condition. More structure is again seen offshore in GBR4 due to its higher resolution there. This structure is even more evident in the temperature solution, however, that aside, the distributions are again generally consistent. The SWR parameterisation in GBR4 was prescribed using a data assimilation parameter estimation technique, which provided an accurate spatially variable parameter field. In the National Model constant default values are currently used, and this is an area that will require parameter optimisation, preferably using similar techniques to eReefs.

Figure 3.3. Sea level and tracers for the national model (top) and eReefs GBR4 (bottom) for 12 Mar. 2017. Note the different times due to time-zones used. Left column is sea level, middle is salinity and right is temperature.

Although resolution is less offshore in the National Model, leading to less meso and sub-mesoscale structure, the opposite is true at the coast and over the reef. Fig. 3.4 shows surface temperature and currents over the reef in an area off the Whitsundays (Fig. 3.4a). GBR4 shows a cooling of water over the reef at this time, and a generally unimpeded southward flow over the reef. The National Model also shows this cooling, with sharper gradients and more structure, but also the flow over the reef shows more variability and is steered by the topography to a greater extent. This is due to more realistic bathymetry resolving of reef flats and channels in the National Model; a result of the 1 km resolution in this region.

Figure 3.4. Surface flow and temperature over the reef off the Whitsundays on 28 Feb. 207. (a) Bathymetry, (b) eReefs and (c) national model.

When persistent south-east winds are present along the Bonney Coast in South Australia, surface offshore Ekman transport drives upwelling at the coast, where the thermocline breaches the surface to result in cold water with an associated strong longshore flow. This is observed in early March 2017 (Fig 3.5a), and is successfully replicated by the model (Fig 3.5b). Upwelling extent and magnitude are comparable to observation, as is the SST within the South Australian gulfs.

Figure 3.5. Surface temperature and currents in the Bonney Coast region demonstrating upwelling (a) IMOS OceanCurrent image and (b) National Model.

The National Model generally replicates the major features of the GBR4 snapshot and dynamics of SE upwelling, and appears to be behaving adequately in its pilot form. Furthermore, some benefits of the model, e.g. lack of cross-shelf boundaries and increased resolution inshore, are apparent. However, the need for further parameter calibration to observations is also apparent.

Sediment / carbon chemistry / optics

The Environmental Modelling Suite contains a selection of suspended sediment, biogeochemical and optical processes. As a demonstration of how these will operate in the National Model v2.95 configuration, we have included:

(1) a 4-layer sediment model with two particle types that can be in suspension or deposited in a sediment layer (Mud and Dust);

(2) a steady-state carbon chemistry model calculated from dissolved inorganic carbon (DIC) and total alkalinity state variables;

(3) air-sea fluxes of oxygen and carbon dioxide that impact DIC and dissolved oxygen concentrations;

(4) a 27-band spectrally-resolved optical model forced by Sun and Moon light and which calculates the downwelling irradiance in the water column and the surface remote-sensing reflectance within the MODIS Bands 3 (blue), 6 (green), and 9 (red) for the calculation of simulated true colour.

Next steps potentially involve adding the pelagic and benthic cycling of the major nutrients through phytoplankton, zooplankton, and benthic communities, with the end goal of creating a national configuration similar to the marine models in the eReefs Project.

4.    Tiling

Once a national mesh is operating, the project aims to partition the mesh into regional tiles that can be operated autonomously while nested with the global model, or ‘stitched’ together via 2-way nesting into a national framework. The partitioning of a national mesh can therefore be achieved such that there is exact overlap of cells across the 2-way interface. This is in essence performing 2-way nesting with no refinement (change in grid size) of the mesh at the nesting interface, and would therefore minimize the spurious short-wave features that are generated and must be controlled when a refinement is present. The technical aspects of tiling the national model therefore revolve around implementation of a robust 2-way nesting system. The approach for this system has been prototyped and assessed (Herzfeld and Rizwi, 2019) and is mature technology that has proven to be suitable for such an application. Decomposition into tiles allows each tile to be run in isolation faster than the national mesh, to result in greater throughput and ultimately a better constrained model. Also, the model parameters can be tuned for each individual tile rather than applying single parameterizations for the whole national domain. An example a national regionalization is shown in Fig. 4.1, with the associated south east Tasmania tile (region 0) shown in Fig. 4.2.

Figure 4.1. Regionalization for tiling.

Figure 4.2. Regional Tasmanian (region 0) tile.

The approach of Herzfeld and Rizwi (2019) used open boundary conditions for information transfer between parent and child models, and found that the transmissive properties of boundary conditions could mitigate specification error issues and allow models to be coupled at the baroclinic (3D) timestep only. This means that sea level is held constant on the boundaries over the 2D time-step, and 2D velocities are the depth average of 3D velocities on the boundary; essentially an average of 2D velocities over the previous 2D time-step. We find that in the proposed application for the National Model, where the separate model tiles are coupled tightly at 1:1 ratios in space and time, that coupling at the barotropic (2D) level is necessary. This is not technically an issue to implement, but does mean that more transfers of information are required (elevation and 2D velocities every 2D time-step). A trial 2-way coupling was performed between tile0 (south-east Australia) and tile5 (southern Australia), with solutions presented in Fig. 4.3. The tile5 solution is shown in Fig. 4.3a, which is coupled to the tile0 solution of Fig 4.3b. The continuous National Model solution is shown in Fig. 4.3c. The parameterization is slightly different between the coupled and continuous solutions, with the horizontal viscosity smaller in the tile0 coupled solution. Nevertheless, solutions are comparable, and again no discontinuities are seen across the interface in the coupled solution. The coupled model also appeared very stable, with no indications of developing instabilities across the interface. The 2-way coupling appears a viable methodology to decompose the national domain into smaller sub-models. While the methodology appears sound, a larger effort will be required to couple all the tiles of the National Model.

Figure 4.3. Two-way nesting of (a) tile5 with (b) tile0. The continuous solution is shown in (c).

5.    Next Steps

The technical barriers to developing a national model are on track, and it is anticipated a scalable pilot model and associated regional tiles will be available by the project’s end. However, the current model achieves maximum average resolution of 2 km at the coast, whereas the target resolution is 1 km. Without efficient parallelization, the development of a 1 km model is unrealistic, therefore a slightly coarser version was developed for addressing the technical issues that were required to be overcome. For the next phase of the project we anticipate to revisit the mesh and produce a highly optimized model at the target resolution.

The activities requiring attention for v3 development are:

1. Mesh development. A higher resolution mesh is to built using existing tools. This is now a routine procedure, however, coastline optimization should be revisited to adequately resolve, or omit, all estuaries, channels and bays around the Australian coastline. Enhanced resolution in certain areas is also to be implemented; initially this is anticipated to be the same as v2 with the addition of waters adjacent to capital cities.
2. The coastline optimization is currently a manual process using rudimentary tools. We require a graphical interface to streamline this process, and some visualization effort is required to be dedicated to developing these tools for the efficient coastline generation required by mesh generation. General progress on MoVE and ParaView will also continue.
3. Once a higher resolution mesh is developed, model stability analysis and validation can commence. A near real-time implementation of v3 will follow, with ongoing parameterization improvement via calibration to observations.
4. The regional tiling can be subject to refinement, and a series of regional tiles produced. These may be used in applications under parallel existing projects, e.g. the GBR tile within eReefs.
5. The parallel processing development will continue, where identification and rectification of bottlenecks will be undertaken. Fine grained and I/O parallelization will be implemented within the code where possible. Algorithms that supplement the dynamic core (e.g. the transport model, 2-way nesting) will also be subject to parallelization.

6.    References

Beaman, R. 2010. Project 3DGBR: a high-resolution depth model for the Great Barrier Reef and Coral Sea. Final Report to Marine and Tropical Sciences Research Facility Final Report, June 2010

Burchard, H., K., O. Peterson and T.P. Rippeth (1998) Comparing the performance of the Mellor Yamada and the k-e turbulence models. J. Geophys. Res., 103, 10,543 – 10,554.

Engwirda, D. (2017) JIGSAW(GEO) 1.0: locally-orthogonal staggered unstructured grid generation for general circulation modelling on the sphere, Geosci. Model Dev. 10 (6) 2117.

Egbert, G.D. and Erofeeva, S.Y. (2002) Efficient inverse modelling of barotropic ocean tides. Journal of Atmospheric and Oceanic Technology, 19(2), 183-204.

Herzfeld, M., Andrewartha, J.R. (2012) A simple, stable and accurate Dirichlet open boundary condition for ocean model downscaling. Ocean Modelling, 43-44, 1-21.

Herzfeld, M., Engwirda, D., Rizwi, F. (2020) A coastal unstructured model using Voronoi meshes and C-grid staggering. Ocean Modelling, 148, 101599. https://doi.org/10.1016/j.ocemod.2020.101599

Kondo, J. (1975) Air-sea bulk transfer coefficients in diabatic conditions. Boundary-Layer Meteorology, 9, 91-112.

Large, W.G., and S. Pond  (1981) Open ocean momentum flux measurements in moderate to strong winds, J. Phys. Oceanogr., 11, 324-336.

Oke, P. R.; G. B. Brassington, D. A. Griffin, A. Schiller, 2008: The Bluelink ocean data assimilation system (BODAS), Ocean Modelling, 21(1-2), 46-70.

Puri, K., Dietachmayer, G., Steinle, P., Dix, M., Rikus, L., Logan, L., Naughton, M., Tingwell, C., Xiao, Y., Barras, V., Bermous, I., Bowen, R., Deschamps, L., Franklin, C., Fraser, J., Glowacki, T., Harris, B., Lee, J., Le, T., Roff, G., Sulaiman, A., Sims, H., Sun, X., Sun, Z., Zhu, H., Chattopadhyay, M., and Engel, C. (2012) Implementation of the initial ACCESS Numerical Weather Prediction system, Aust. Met. Oc. Journal, 63, 265-284.

Ringler, T., Petersen, M., Higdon, R. L., Jacobsen, D., Jones, P. W., and Maltrud, M. (2013) A multi-resolution approach to global ocean modeling. Ocean Modelling, 69(C), 211–232. doi:10.1016/j.ocemod.2013.04.010

Ringler, T., Petersen, M., Higdon, R. L., Jacobsen, D., Jones, P. W., and Maltrud, M. (2013) A multi-resolution approach to global ocean modeling. Ocean Modelling, 69(C), 211–232. doi:10.1016/j.ocemod.2013.04.010

Thuburn, J. (1995) Multidimensional flux-limited advection schemes. J. Comput. Phys., 123, 74-83.

Thuburn, J., T. Ringler, J. Klemp and W. Skamarock, 2009: Numerical representation of geostrophic modes on arbitrarily structured C-grids, Journal of Computational Physics, 2009: 228 (22), 8321-8335. doi:10.1016/j.jcp.2009.08.006

Zillman, J.W. (1972) A study of some aspects of the radiation and heat budgets of the southern hemisphere oceans. Bureua of meteorology, Meteorological study no. 26, Australian Govt. Pub. Service, Canberra.

7.    Downloads

Progress report 2020 Download

Stakeholder ppt