MTSRF (GBR, 2013-2016)

MTSRF Climate Downscaling Project

This project is conducted within the Marine and Tropical Sciences Research facility (MTSRF);


Coral bleaching is now perceived as a major threat for coral reef ecosystems, especially since the major ENSO-related 1998 event (Hoegh-Guldberg(1999), Wilkinson 2000). Some of the impacted reefs have not recovered, resulting in negative changes in fish and benthic community structures, biodiversity, productivity and fisheries resources (e.g. Chabanet 2002).

Predicting the vulnerability to bleaching at the scale of few tens of kilometers is of importance when designing monitoring programs or to design Marine Protect Areas (MPA) since resistance to bleaching stress is a desirable property of areas intended to protect biodiversity.

Coral bleaching is mostly driven by positive anomalies in Sea Surface Temperature (SST) (Wellington et al. 2001). Consequently, SST monitored from space using the AVHRR sensors have been used for nearly a decade as relatively efficient indicators of potential bleaching (Liu et al. 2003). However, ultimately, it is the local and regional hydrodynamic regime that is the primary controlling factor since SST is driven by mixing conditions in the upper water column (Brown 1997). Topography, wind, low frequency currents and tide regimes are critical information to understand bleaching patterns.

The GBR is a 2000 km long reef and lagoon complex. It is a UNESCO World Heritage Area covering 347,800 sq. km, sandwiched between the Australia’s northeastern coastline and the Coral Sea. It includes an estimated 2800 fringing, patch and shelf reefs organized in a matrix filled by an immense lagoon of average depth ~35m. It supports $AUS250 million of commercial, ships and tourism industries (Wachenfield et al. 1998). Heron Island on the southern GBR is a unique experimental site because long term biological obervations due to the presence of the research station and continuous physical oceanographic observations since 2004 allow this project to leverage off this archive of information.

A pilot nested modelling suite has been developed to downscale global scale model simulations to the reef scale, for the purpose of investigating climate change on coral bleaching. The focus for coral bleaching analysis was chosen to be the vicinity of Heron Island, hence a nesting strategy was required that enabled major oceanographic processes in the GBR to be propagated into a high resolution model domain surrounding this site. The global model simulations comprise three climate scenarios corresponding to the climate in 2070 if carbon emissions are:

– Capped at today’s sea levels,
– Extrapolated from current trends,

Using the nested model suite to downscale climate scenarios to the reef scale (< 100 m resolution), the effect of these future climates on coral bleaching may be investigated.

Regional Oceanography

The general circulation of the Great Barrier Reef region is described by Burrage et al (1996). The South Equatorial Current flows into the Coral Sea where it splits into two as it approaches the northeast Australian shelf. A southern arm flows south, trapped against the slope, to form the East Australian Current (EAC). The northern arm flows on the shelf edge in the northern Great Barrier Reef (GBR) to form a semi-closed cyclonic eddy in the Coral Sea. In the northern GBR upwelling due to fluctuations in geostrophic currents or northeast monsoon winds transports cool water along the bottom of the central GBR shelf to the seaward edge of the GBR lagoon. This cool water does not create a surface signature. However, the cool water may reach the surface via upwelling due to internal or barotropic shelf edge tides, island wakes, topographically induced eddies or bottom Ekman pumping.
The circulation and thermal characteristics of the model region have been studied by Griffin et al (1987). The Capricorn Channel is associated with large tides of amplitude approaching 4 m at the coast. The tide is of semidiurnal nature in the north, progressing to equal diurnal and semidiurnal behaviour in the south near Lady Elliot Is. Flood tide directions have a strong alongshore component in the GBR lagoon, becoming stronger and later further north. Non-tidal circulation is complex and highly variable. Generally, circulation off the shelf break is dominated by pulses of north-westward flow greater than 0.3 ms-1, having a period of 6 – 10 days. This large flow is due to the cumulative effect of several baroclinic coastally trapped wave modes (Griffin and Middleton, 1986) and gives rise to a mean north to north-westward current predominantly parallel to the isobaths. The origin of these coastally trapped waves is further south beyond Fraser Island. Local wind stress contributes to circulation on the shelf, and to a lesser extent to mean flow on the shelf break. Griffin et al. (1987) postulated the presence of a large cyclonic eddy located at the southern mouth of the Capricorn Channel, generated by the EAC. This results in a north-westward backflow on the shelf-break in the vicinity of Lady Musgrove Is. The appearance and disappearance of this eddy in response to EAC fluctuations may give rise to a long period (90 day) variation in the mean current on the shelf-break.

Near the coast the temperature does not fluctuate in response to the tide. In shallow water near reef edges diurnal fluctuations of up to 4oC were observed by Griffin et al (1987) and attributed to solar heating effects. Further offshore a tidal oscillation of up to 2oC was observed on a flood spring tide. Larger temperature oscillations are observed in deeper waters and are attributed to tidal pumping. The temperature lags the offshore current component by 90o (temperature minima occur at the end of the incoming tide). Griffin et al. (1987) show that a tidal current of 0.25 ms-1 may induce a vertical displacement of 60m resulting in ~1.5oC temperature variation. The thermocline also oscillates on longer timescales (6 – 10 days) on the shelf break and slope. These fluctuations were shown by Griffin and Middleton (1986) to be related to the alongshore current due to coastally trapped waves resembling an internal Kelvin wave.

The Model

The model used in this study, SHOC (Sparse Hydrodynamic Ocean Code; Herzfeld et al., 2006b, Herzfeld and Waring, 2006) is a general purpose model applicable to scales ranging from estuaries to regional ocean domains, and has been successfully applied to a variety of applications encompassing these scales (e.g. Walker (1997), Herzfeld et al. (2006a)). This model is a three-dimensional finite difference hydrodynamic model based on the primitive equations. Outputs from the model include three-dimensional distributions of velocity, temperature, salinity, density, passive tracers, mixing coefficients and sea level. Inputs required by the model include forcing due to wind, atmospheric pressure gradients, surface heat and water fluxes and open boundary conditions (e.g. tides). A schematic of the forcing mechanisms used in SHOC is presented in Figure 3.1. SHOC is based on the three dimensional equations of momentum, continuity and conservation of heat and salt, employing the hydrostatic and Boussinesq assumptions. The equations of motion are discretized on a finite difference stencil corresponding to the Arakawa C grid.
The model uses a curvilinear orthogonal grid in the horizontal and a choice of fixed ‘z’ coordinates or terrain-following s coordinates in the vertical. The ‘z’ vertical system allows for wetting and drying of surface cells, useful for modelling regions such as tidal flats where large areas are periodically dry. The current implementation of the model uses z-coordinates. The bottom topography is represented using partial cells. SHOC has a free surface and uses mode splitting to separate the two dimensional (2D) mode from the three dimensional (3D) mode. This allows fast moving gravity waves to be solved independently from the slower moving internal waves allowing the 2D and 3D modes to operate on different time-steps, resulting in a considerable contribution to computational efficiency. The model uses explicit time-stepping throughout except for the vertical diffusion scheme which is implicit. A Laplacian diffusion scheme is employed in the horizontal on geopotential surfaces. Smagorinsky mixing coefficients may be utilized in the horizontal.

Schematic of forcing mechanisms in SHOC.
Figure 3.1: Schematic of forcing mechanisms in SHOC.

The ocean model can invoke several turbulence closure schemes, including k-e, Mellor-Yamada 2.0 and Csanady type parameterizations. A variety of advection schemes may be used on tracers and 1st or 2nd order can be used for momentum. The model also contains a suite of open boundary conditions, including radiation, extrapolation, sponge and direct data-forcing. A generous suite of diagnostics in included in the model.

The equations of motion SHOC is based upon are similar to those described in Blumberg and Herring (1987). Details of the model equations can be found in Herzfeld (2006a).

Nested Modelling Strategy

It is desirable that all the dominant physical processes are captured in a local model domain of Heron Island. For this to be accomplished a three grid nesting strategy was employed (Figure 4.1). The large scale (regional) grid is designed to capture EAC intrusions onto the shelf and resolve the northward propagation of coastally trapped waves south off Fraser Island. These internal waves are hypothesised to collide with Fraser Island and degenerate into higher frequency modes that continue to propagate northwards. The resolution of global model products (e.g. BRAN reanalysis) is too coarse (10 km) to adequately resolve these higher frequency waves, hence the resolution of the regional model (3 – 6 km in this region) must fulfil this role. This grid has a minimum depth of 2 m imposed at the coast, and maximum depth is 2500 m. An intermediate model of the Capricorn Group acts as a nesting vehicle to allow boundary conditions to be downscaled to a resolution applicable to driving the local Heron Island grid, and incorporates the influence of tides. Resolution of this grid varies from 800m to 4km, with a maximum depth in the domain of 320 m. Finally the local grid is a cyclic polar grid that maintains high resolution (down to 12m) around Heron Island and increasing to ~2km at the grid edge to accept boundary forcing from the intermediate grid. The maximum depth of this grid is 44 m.The maximum depth of the models impacts on the gravity wave speed, hence the time step used for the 2-D mode and consequently the run time of the model. The regional model has 3D/2D time-steps of 200/4 s respectively. The grid size if 100 x 270 with 38 layers in the vertical, with ~33% of the grid containing wet cells. This results in a run-time ratio of ~110:1 using parallel processing over two partitions. The intermediate grid is considerably slower. The size of this grid is 110 x 170 with 34 vertical layers, with ~68% wet cells. The 3D/2D time-steps are 104/4 s, resulting in run-time ratios of ~40:1 when parallel processed over 4 partitions. Finally the local grid has a size of 70 x 37 with 22 vertical layers (~42% wet cells) and 3D/2D time-steps of 6/1 s. Without parallel processing this grid has a run-time ratio of ~25:1, making it suitable only for shorter term seasonal simulations.


All grids use ‘z’ vertical discretization with exponentially increasing grid spacing near the surface and constant spacing at depth. Surface layer thickness is 1, 0.3 and 0.2 m for the regional, intermediate and local grids respectively. The bathymetry for all grids is smoothed with a 9 point convolution filter, and a maximum gradient of 0.05 is imposed. All models use the higher order upwind advection scheme if Van Leer (1979).

The regional model is forced with global model products. Alternatively, altimeter / climatology based products may be used (e.g. synTS) but since these products rely on climatology to project surface values to depth they are not able to represent internal oscillations of the thermocline. The summer period Dec 05 to Mar 06 will be the focus of coral bleaching analysis, hence of model calibration. During this period a field program was in place by AIMS to produce a sufficient quantity of measurements for comparison to model output. The regional model was nested in BRAN2.1 for this period; these outputs considered the best global product to date suitable for forcing the nested suite. Atmospheric forcing products (wind, pressure, heatflux) are supplied by the MesoLAPS ( atmospheric model. This modelled product is preferable to spatially interpolated meteorological data supplied by BOM weather stations for surface forcing, since its spatial detail is superior. MesoLAPS provides wind, mean sea level pressure, short wave radiation, air temperature and dew point temperature. From these variables bulk schemes (Kondo, 1975) are used to compute sensible and latent heat fluxes, and black body radiation may be used to compute long wave radiation (Zillman, 1972). The sum of these and short wave input provides a net heat budget. The heat flux components and net heat budget are displayed in Figure 4.2.

A tide model was prepared by Bode et al. (1997) which provided amplitudes and phases of 19 tidal constituents (Table 4.1) at a resolution of about 2 km. This tidal model was not large enough to cover the regional domain, hence the regional model was simulated without tides, and the tidal signal was superimposed on the intermediate model’s open boundaries in conjunction with sea level derived from the regional model. It is observed from Table 3.1 that the dominant constituents in the region are those due to M2, S2 ,K1, O1 and N2. The tide is mixed mainly semi-diurnal in character, which may be quantified by the form factor F = ratio of diurnal to semi-diurnal amplitudes (F = K1+O1 / M2+ S2). In the case of Heron Island F~0.36, verifying that the tide falls into the semi-diurnal mixed category (0.25 < F < 1.5).

Model nesting strategy.
Figure 4.1: Model nesting strategy.


Heat budget at heron island from MesoLAPS model.
Figure 4.2: Heat budget at heron island from MesoLAPS model.


Tidal harmonics for Heron Island.
Table 4.1: Tidal harmonics for Heron Island.

Field Measurements
Data from the oceanographic array in the Capricorn Bunker Group is to be used to calibrate and force the model. The array over the summer of 2005-2006 had a total of sixty-nine instruments deployed for this experiment throughout the Heron Island region. The deployments and recoveries took place during four field trip periods between the months of August 2004 and April 2006. Many instruments were positioned in two or more locations during the study. The following figures show the locations of the main array in the vicinity of Heron Island.

The array consisted of current meter moorings designated M0, M1, M2, M3, M4 and MB (Figure 5.1). MB, north of Wistari Reef consisted of a surface meteorological buoy that also sampled incoming and outgoing short and long wave radiation fluxes. All the other moorings consisted of current meter profilers and temperature loggers at various depths throughout the water column. In addition to the moorings, temperature transect were deployed around Heron Island (Figure 5.2). These were designated T7, T9, T10 and T11. T9 effectively comprised a transect across the eastern end of the reef from south to north.

Cross shelf current meter and temperature transect.
Figure 5.1: Cross shelf current meter and temperature transect.


Heron island temperature transect locations.
Figure 5.2: Heron island temperature transect locations.


Magnification of the boxed area at the western end of Fig. 5.2. This is the channel leading to the Heron island pier and shows the location of transect T7 on the side of reef and a current meter in the Heron island Channel (HIC).
Figure 5.3: Magnification of the boxed area at the western end of Fig. 5.2. This is the channel leading to the Heron island pier and shows the location of transect T7 on the side of reef and a current meter in the Heron island Channel (HIC).


Fluxes, temperature, wind stress and sea level from the meteorological buoy.
Figure 5.4: Fluxes, temperature, wind stress and sea level from the meteorological buoy.


Sea level and temperature transect data from Heron Reef Slope in the Wistari Channel, T7.
Figure 5.5: Sea level and temperature transect data from Heron Reef Slope in the Wistari Channel, T7.

Model Output
Models representing the GBR-scale, Capricorn Bunker Group and reef-scale around HeronIsland models were re-configured to simulate this period using the described forcing data to produce stable simulations. Examples of output for all nesting scales are displayed in Figures 6.1 – 6.3.

An initial comparison between the model and in situ observations over a 6 week period are presented in Figures 6.4 and 6.5. Figure 6.4 shows the temperature traces of water at the bottom of the T7 transect, approximately 17m deep. The model correlates quite well with the long term trend of the gradual summer warming although the data precedes the model with a minor cooling period. The model does eventually cool and continues to track the observations well.

The delay in the cooling episode is perhaps due to a combination of the relatively coarse forcing of the modelled large scale circulation (updated daily) and the atmospheric forcing data. Nevertheless the main long term warming trends are quite well represented in this topographically complex location.

Of note is the significant higher frequency tidal signal evident in the observations. These have significantly higher amplitudes than what is found in the model.

Figure 6.5 shows the equivalent time series of sea level (n.b. the data is offset due to the different datums). This reveals the tidal forcing in the model is under-representing the true strength of the tide in this area. It is expected that once the tidal forcing is made more realistic, the modelled daily fluctuations in temperature will be much more closely correlated with the observations. Comparisons of modelled and measured velocity are shown in Figure 6.6

There needs to be some further calibration of the model in order to improve the comparisons. This will be achieved by improving the forcing parameters by including improved observational data, comprising actual observed air-sea energy fluxes and further tidal forcing improvements. The tides are responsible for the observed significant daily variability in temperature and so this is a critical component.

The first downscale scenario has been propogated through all nesting levels. The surface temperature and currents for the Heron island model are shown in Figure 6.7, from which it is seen that water is significantly warmer than that displayed in Fig. 6.3.

Surface temperature and currents derived from the regional model.
Figure 6.1: Surface temperature and currents derived from the regional model.


Surface temperature and currents derived from the Capricorn model.
Figure 6.2: Surface temperature and currents derived from the Capricorn model.


Surface temperature and currents derived from the Heron model.
Figure 6.3: Surface temperature and currents derived from the Heron model.


Bottom temperature comparison between model (blue) and observed (red) in Wistari Channel location T7.
Figure 6.4: Bottom temperature comparison between model (blue) and observed (red) in Wistari Channel location T7.


Sea level comparison between model (blue) and observed (red) at Wistari Channel location T7.
Figure 6.5: Sea level comparison between model (blue) and observed (red) at Wistari Channel location T7.


Velocity comparison between model (blue) and observed (red) at Wistari Channel location T7.
Figure 6.6: Velocity comparison between model (blue) and observed (red) at Wistari Channel location T7.
Surface temperature and currents around heron Island from the 1st downscaled scenario.
Figure 6.7: Surface temperature and currents around heron Island from the 1st downscaled scenario.



Berkelmans R, Oliver J (1999) Large scale bleaching of corals on the Great Barrier Reef. Coral Reefs 18:55-60
Blumberg, A.F. and J. Herring (1987) Circulation modelling using orthogonal curvilinear coordinates, in Three-Dimensional Models of marine and Estuarine Dynamics, Ed. J.C.J. Nihoul and B.M. Jamart, Elsevier.
Bode, L., Mason, L.B. & Middleton J.H. (1997) Reef parameterisation schemes with applications to tidal modelling. Prog. Oceanogr. 40, 285-324.
Brown BE (1997) Adaptations of reef corals to physical environmental stress. Advances in Marine Biology 31:221-299
Burrage, D.M., C. R. Steinberg, W. J. Skirving and J. A. Kleypas (1996)Mesoscale circulation features of the Great Barrier Reef lagoon inferred from NOAA satellite imagery. Remote Sens. Environ., 56, 21 – 41.
Chabanet P (2002) Coral reef fish communities of Mayotte (western Indian Ocean) two years afer the impact of the 1998 bleaching event. Marine Freshwater Research 53:107-113
Done T, Turak E, Wakeford M, De’ath G, Kininmonth S, Wooldridge S, Berkelmans R, Van Oppen M (2003) Testing bleaching resistance hypotheses for the 2002 Great Barrier Reef bleaching event. Unpublished report to The Nature Conservancy. Australian Institute of Marine Science: 106 p.
Griffin, D. A. and J. H. Middleton (1986) Coastal trapped waves behind a large continental shelf island, southern Great Barrier Reef. J. Phys. Oceanogr., 16, 1651 – 1664.
Griffin, D. A., J. H. Middleton and L. Bode (1987) The tidal and longer-period circulation of Capricornia, souther Great Barrier Reef. Aust. J. Mar. Freshw. Res., 38, 461 – 474.
Herzfeld, M. (2006) An alternative coordinate system for solving finite difference ocean models. Ocean Modelling, 14, 174 – 196.
Herzfeld, M., J. Andrewartha, P. Sakov and I.T. Webster (2006a) Numerical hydrodynamic modeling of the Fitzroy Estuary. CRC for Coastal Zone, Estuary and Waterway Management Technical report 38, 95pp.
Herzfeld, M., J. Waring, J. Parslow, N. Margvelashvili, P. Sakov and J. Andrewartha (2006b) SHOC, sparse hydrodynamic ocean code, Scientific manual. CSIRO Marine Research, 120pp.
Herzfeld, M., and J. R. Waring (2006) SHOC, sparse hydrodynamic ocean code, User manual, CSIRO Marine Research, 118pp.
Hoegh-Guldberg O (1999) Climate change, coral bleaching and the future of the world’s coral reefs. Marine and Freshwater Research 50:839-866
Kantha, L.H. and C.A. Clayson (2000) Small scale processes in geophysical fluid flows. Vol. 67 of International Geophysics Series, Academic Press, 2000.
Kondo, J. (1975) Air-sea bulk transfer coefficients in diabatic conditions. Boundary-Layer Meteorology, 9, 91-112.
Liu G, Strong A, Skirving W (2003) Remote sensing of Sea Surface Temperature during the 2002 (Great) Barrier Reef coral bleaching. EOS Trans. AGU 84:137, 144
Mann K.H. and J.R.N. Lazier (1996) Dynamics of Marine Ecosystems : Biological-Physical interactions in the ocean. 2nd Ed. Blackwell Science.
Steinberg C (2007) Impacts of climate change on the physical oceanography of the Great Barrier Reef. Chapter 3, Climate Change and the Great Barrier Reef, eds. Johnson JE and Marshall PA. Great Barrier Reef Marine Park Authority & Australian Greenhouse Office, pp 51-74.
Steinberg, C., Mahoney, M., Andréfouët. S., and Skirving, W. (2002)Recognising cool reefs in ocean hotspots: remote sensing and modeling approaches. Transforming Coral Reef Conservation in the 21st Century , Coral Beaching Targeted Research Working Group. The Nature Conservancy and Conservation International 18-19 March 2002, Townsville, Australia.
Wachenfield D (1998) State of the Great Barrier Reef World Heritage Area. 140 p. GBRMPA.
Walker, S.J. (1997) Hydrodynamic models of Port Phillip Bay. CSIRO Port Phillip Bay Environment Study Technical Report no. 38, Melbourne.
Wellington GW, Glynn PW, Strong AE, Navarrete SA, Wieters E, Hubbard D (2001) Crisis on coral reefs linked to climate change. EOS Trans. AGU 82:1,5
Wilkinson C (2000) Status of the coral reefs of the world: 2000. Australian Institute Marine Sciences, 363 pp
Van Leer, B. (1979) Towards the ultimate conservative difference scheme. V: a second order sequel to Godanov’s method. J. Comput. Phys., 32, 101-136.
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.