Assessment of the Accuracy of Recent Empirical and Assimilated Tidal Models for the Great Barrier Reef, Australia, Using Satellite and Coastal Data

Abstract

The latest satellite and in situ data are a fundamental source for tidal model evaluations. In this work, the satellite missions TOPEX/Poseidon, Jason-1, Jason-2 and Sentinel-3A, together with tide gauge data, were used to investigate the performance of recent regional and global tidal models over the Great Barrier Reef, Australia. Ten models, namely, TPXO8, TPXO9, EOT11a, HAMTIDE, FES2012, FES2014, OSUNA, OSU12, GOT 4.10 and DTU10, were considered. The accuracy of eight major tidal constituents (i.e., K₁, O₁, P₁, Q₁, M₂, S₂, N₂ and K₂) and one shallow water constituent (M₄) were assessed based on the analysis of sea-level observations from coastal tide gauges and altimetry data (TOPEX series). The outcome was compared for four different subregions, namely, the coastline, coastal, shelf and deep ocean zones. Sea-level anomaly data from the Sentinel-3A mission were corrected using the tidal heights predicted by each model. The root mean square values of the sea level anomalies were then compared. According to the results, FES2012 compares more favorably to other models with root mean square (RMS) values of 10.9 cm and 7.7 cm over the coastal and shelf zones, respectively. In the deeper sections, the FES2014 model compares favorably at 7.5 cm. In addition, the impact of sudden fluctuations in bottom topography on model performances suggest that a combination of bathymetric variations and proximity to the coast or islands contributes to tidal height prediction accuracies of the models.

1. Introduction

The Great Barrier Reef (GBR), Australia, is the largest reef system in the world. It comprises over 2500 coral reefs and stretches for more than 2000 km along the northeastern Australian coastline, ranging in width from 100 to 200 km [1]. The area has an average depth of 35 m in its inshore waters, and, on the outer reefs, the continental slope depths extend to over 2000 m, but it is categorized overall as a shallow water region. In comparison to non-reef zones, reefs can exert noticeable effects on tides and other large-scale flows in certain parts of the GBR, leading to complex dynamic processes in the region [2]. Many investigations of the dynamic phenomena, sea level variations, tides and surface currents in the GBR have been conducted [2,3,4,5,6,7,8,9,10,11,12].

The resonant and near-resonant responses of tidal waves and the significance of the nonlinearity in dynamic equations increases the complexity of tidal analysis in coastal regions [13,14]. In the GBR, several factors account for the complex tidal behavior. The bottom drag coefficient, which is relevant to major tidal components, varies noticeably over reef zones [15]. There are also considerable changes in bottom topography [16]. The variable bathymetry and topography of this region, as well as the existence of reefs, renders it difficult to capture a realistic image of the hydrodynamics of the area, including the tides [17]. The wavelength of the major tidal constituents is on the order of 1000 km at typical reef depths. Consequently, individual reefs in this region may not exert a significant impact on tidal responses [2]. The tidal response is more likely to be affected by a long chain of reefs with a very large variation in bottom topography in the southern and central GBR [4]. For example, the dense coral reefs directly offshore from Broad Sound hinder the normal flow of tidal currents, for which the high tide at Broad Sound consists of two currents flowing from the north and from the south. This feature magnifies the amplitude of the semidiurnal constituents [6]. In the central GBR, a study by Andrews and Bode [7] showed that the two major semidiurnal constants, M₂ and S₂, are amplified by 50%. The surface and tidal currents are trapped in gyres, waves and eddies, and they are split into small branches mainly due to the bottom topographic variations of the reefs [4,12].

In addition, the lack of sufficient information regarding the shallow water constituents (including both compound and overtide components) in existing tidal models leads to the absence of a comprehensive picture of tides and circulations in the area [18]. These tidal components, which are the consequences of interactions between bottom topography and major tidal constants, exhibit large variations in the presence of the reefs. For example, constituents M₄ and MS₄ are higher-frequency harmonics with amplitudes that are usually much smaller than M₂. However, in locations where shallow water exerts impacts on tides, their amplitudes are comparable with the M₂ constituent [13]. This distortion due to shallow water tidal components can increase when the drag coefficient, and subsequently the bottom friction, changes substantially due to the coral reefs and changes in bathymetry [15]. Accordingly, to comprehensively present the tides and circulations in the GBR, it is essential to obtain sufficient information about the compound and overtide shallow water components.

Global data from satellite altimetry and from regional tidal models have seen great improvement. However, these data are known to be less reliable in coastal zones than in the open ocean because of inaccurate geophysical corrections and noisier radar impulse responses [19]. The development in recent years of more efficient radar retracking algorithms and improved geophysical models has led to more accurate sea surface height (SSH) over the coastal zones. Previous studies evaluating recent tidal models (both regional and global) have found that the recent models are more accurate than their predecessors [20,21,22,23]. However, the performance of these models in shallow water and coastal zones is inferior to their performance in the open ocean in terms of the accuracy of tidal constituent estimations and, consequently, the accuracies of tidal height predictions [24]. Given that tidal phenomena in the GBR are typically more complex than those in other coastal zones, detailed model assessments in this region are a precursor to further tidal studies by revealing the strengths and weaknesses of existing models.

In this study, ten recent empirical, semi-empirical and assimilation tidal models were included for assessment [21]. The first assessments reflected the ability of models to estimate tidal constants and the contribution of other aspects of model development such as spatial resolution. In the second assessment, the models were compared with Sentinel-3A data. Sentinel-3A data were not assimilated in any of the aforementioned models and therefore provide criteria for an independent comparison. Section 2 briefly describes the study area, data and the tidal models used in this study. In Section 3 and Section 4, the evaluation of the models for tidal constant estimations and tidal height predictions is shown. Section 5 presents a discussion regarding the performance of the models and subsequent conclusions.

2. Data and Study Area

The study area is bounded by latitudes 10° S to 29° S and longitudes 142° E to 160° E, covering a total area of ~3,000,000 km², which covers the GBR on the northeastern coast of Australia (Figure 1). In the GBR, coral reefs extend from the coastline to the shelf zone. Figure 1 shows the bathymetry based on the GBR100 model available from https://www.deepreef.org/bathymetry/65-3dgbr-bathy.html. This model utilized the latest data sources, including data from ship-based multibeam, single-beam and echo-sounder surveys; airborne LiDAR bathymetry surveys; and satellite remote sensors [25]. It is currently the most accurate bathymetry model, with a high spatial resolution of 3.5″ × 3.5″ in the study area (Figure 1). In this study, the area of interest was divided into four zones depending on the ocean depth as follows:

• Coastline
• Coastal zone with a depth of 0–40 m
• Shelf zone with a depth of 40–400 m
• Deep ocean zone with a depth of over 400 m

The coastline consists of only tide gauge stations, while the other zones consist of only altimetry along-track locations.

2.1. Tide Gauge Records

Sea-level observations recorded at 18 tide gauges between 1991 and 2015 (

Table 1. The geographic positions and observation record periods of the coastal tide gauges.

Station No.	Tide Gauge Name	Latitude (Degrees)	Longitude (Degrees)	Time Duration
				Start Date	End Date
1	Bowen Storm	−20.01	148.25	1 January 1996	31 December 2015
2	Bundaberg	−24.76	152.38	1 January 1996	31 December 2015
3	Cairns Storm	−16.93	145.78	1 January 1996	31 December 2010
4	Cape Ferguson	−19.26	147.05	1 January 1996	31 December 2015
5	Cardwell	−18.25	146.03	1 September 1991	31 December 2015
6	Clump Point	−17.83	146.1	1 January 1997	31 December 2015
7	Cook Town	−15.45	145.25	1 January 1996	31 December 2015
8	Fishermans Landing	−23.78	151.16	1 January 2002	31 December 2015
9	Gladstone	−23.83	151.25	1 January 1996	31 December 2015
10	Mackay	−21.1	149.21	1 January 1996	31 December 2006
11	Mourilyan	−17.6	146.11	1 January 1996	31 December 2015
12	Port Alma	−23.58	150.86	1 January 1996	31 December 2015
13	Rosslyn Bay	−23.16	150.8	1 January 2001	31 December 2015
14	Shute Harbour	−20.28	148.78	1 January 1996	31 December 2015
15	South Trees	−23.85	151.3	1 January 1996	31 December 2015
16	Thursday Island	−10.56	142.21	1 January 1996	31 December 2015
17	Townsville	−19.25	146.83	1 January 1996	31 December 2011
18	Urangan	−25.3	152.91	1 January 1996	31 December 2015

) were used. The data have a high accuracy of 1 cm under all weather conditions [26], and their sample frequencies are either 1 h or 0.5 h. They can be used to extract all tidal components including both astronomical and compound/over tides due to their high temporal resolution. The tide gauge data are available from the Australian Bureau of Meteorology (BoM, https://www.bom.gov.au/oceanography/projects/abslmp/data/) and from Maritime Safety Queensland (https://www.msq.qld.gov.au/Tides/Open-data). Table 1 shows that Mackay Station has been operating for the minimum time length of 10 years, while the other tide gauges have had up to 19 years of sea level observations.

2.2. Altimetry Data

Figure 2 shows TOPEX/POSEIDON (T/P) and Sentinel-3A orbits up to 400 km from the coast covering the entire GBR region. The orbit of T/P was specifically designed for extracting ocean tides, and the combined T/P+Jason-1+Jason-2 time series has a relatively short repetition period of approximately 10 days and a typical accuracy of 2 cm [27]. However, the groundtrack resolution is more than 300 km, which is problematic for tidal modeling in the coastal zone.

In this study, delayed time corrected sea surface height (CorSSH) data (1996–2014) provided by Archiving, Validation and Interpretation of Satellite Oceanographic Data (AVISO, https://www.aviso.altimetry.fr) were used. This product was chosen due to its specific editing criteria, which include optimal geophysical corrections in the coastal areas [28]. The CorSSH products provide SSH from multiple altimetry satellites. The geophysical and environmental corrections applied by the data supplier are dry troposphere, wet troposphere, ionosphere, sea state bias, ocean tide and loading tide, solid earth tide, pole tide and dynamic atmospheric correction (DAC) corrections [29]. Since tides are of interest, ocean tide and loading tide corrections are added back to the SSHs. The sea level anomaly that includes both ocean and load tides (SLAt) is calculated as:

$$SLAt = SSH - MSS$$

(1)

where the SSH is from CorSSH, and includes both ocean tides and the loading tide, and MSS is the mean sea surface from the CNES-CLS11 [30].

To evaluate the performance of tidal models in predicting tidal heights, an independent dataset that has not been used in model implementation should be considered. The Sentinel-3A altimetry mission, launched in February 2016 by the European Space Agency (ESA), is a part of the Sentinel missions designed to provide measurements for monitoring the global environment [31]. Sentinel-3A has a repeat period of 27 days and has been designed to observe the sea surface with an accuracy of ~2–3 cm [32]. Sentinel-3A operates in a sun-synchronous orbit and thus fails to detect sun-related constituents, such as S₁ and S₂. However, this is not a major concern for this study, as this mission was used only to assess the ability of models for tidal height predictions. Sentinel-3A SSHs are available from the Radar Altimeter Database System (RADS, rads.tudelft.nl). The software available with RADS [33] was used to extract SLA time series over the GBR region (Figure 2). At the time of this study, 40 cycles of Sentinel-3A SLA, corresponding to three years of data (2016–2019), were available.

2.3. Tidal Models Used

To understand the capabilities of existing barotropic tidal models for tidal constant extraction and prediction, ten models were assessed in this study (

Table 2. The tidal models used in this study.

Model	Type*	Scale#	Tide Gauge	Resolution	Altimeter Data$	Authors
DTU10	S	G	No	1/8°	T/P, J1, J2, EN, E2	Cheng and Andersen (2011)
TPXO8	A	G	Yes	1/30°	T/P, J1, J2, EN, E2	Egbert (2010)
FES2012	A	G	No	1/16°	T/P, J1, J2, EN, E1,E2	Lyard et al. (2006, Updated)
FES2014	A	G	Yes	1/16°	T/P, J1, J2, EN, E1,E2, TPN-J1N	Lyard et al. (2006, Updated)
EOT11a	S	G	No	1/8°	T/P, J1, J2, EN, E2	Savcenko and Bosch (2012)
HAMTIDE	A	G	No	1/8°	T/P, J1, J2, EN, E2	Taguchi et al. (2014)
OSU12	E	G	No	1/4°	T/P,J1,GFO,EN	Fok (2012)
OSUNA	A	R	No	1/24°	T/P, J1, J2	Ray et al. (2011)
GOT4.10	S	G	No	1/2°	J1, J2, GFO, E1, E2, EN, ICEsat	Ray (1999,Updated)
TPXO9	A	G	Yes	1/30°	T/P, J1, J2, EN, E2	Egbert (2010, Updated)

* E, Empirical; S, Semi-empirical; A, Assimilation. ^# G, global; R, regional. ^$ J1, Jason1; J2, Jason2; EN, Envisat; E1, ERS1; E2, ERS2; TPN, Topex Interleaved Orbit; J1N, Jason1 Interleaved Orbit; GFO, GEOSAT Follow-on.

). These models represent the state of recent improvements in empirical and assimilation ocean tidal modeling that have been summarized in [21] and are briefly introduced in this paragraph.

The ocean tide models can be grouped into three categories: empirical models, such as OSU12, which rely only on analyzing SSH observations made by altimetry missions and by tide gauges; semi-empirical models (DTU10 and EOT11), in which a barotropic hydrodynamic model is used to remove tides from the SSH measurements and to then analyze the time series of residuals for tidal constituent corrections; and barotropic hydrodynamic models (OSUNA, TPXO8, TPXO9, FES2012, FES2014 and HAMTIDE) that are constrained using empirical observations through different assimilation approaches.

The models included in this study are briefly described below.

• DTU10 is a global semi-empirical model by Cheng and Andersen [34]. It uses 17 years of observations from the T/P, Jason-1 and Jason-2 missions as well as observations from ERS2 and Envisat. To implement this model, tidal residuals are calculated using FES2004 [35]. Tidal constants are estimated using the response method on tidal residuals. DTU10 benefits from a depth-dependent method for gridding, increasing its ability to estimate tidal constants over areas with highly varying bathymetry as well as on the coastline [21].
• TPXO8 is an assimilation model implemented by [36]. Altimetry and coastal tide gauge data are assimilated into shallow-water equations using a representer method. This model benefits from the GEBCO bathymetry model and from regional bathymetry data and has been developed based on the Oregon State University Tidal Inversion Software (OTIS).
• FES2012 is an extension to the finite-element solution (FES) series started by [37]. It is an assimilation model in which the T/P series, ERS and Envisat data have been assimilated using an improved representer approach. This global model includes 32 tidal constituents.
• FES2014 is the latest version of the FES series and usually shows higher accuracy in shallow water zones because of finer bathymetry models and an optimized assimilation scheme [38]. Moreover, its finite element grid is twice the size that used in FES2012. In addition to altimetry, tide gauge data are also assimilated in this model.
• EOT11a is based on FES2004 with tidal constituent corrections deduced from a harmonic analysis of satellite altimetry data. This model has focused on improvements over shallow water regions by using two steps in the preprocessing phase. First, the most recently updated retracked data are used (over coastal zones), and, second, the same updated geophysical corrections for all the implemented satellite data are applied [39].
• HAMTIDE was developed by Taguchi et al, [40] and is an assimilation global ocean model that uses linearized hydrodynamic equations on a global scale. The constraints have been estimated from EOT11a and assimilated into linearized shallow water equations. In this work, the HAMTIDE11a version of the model was assessed.
• OSU12 is an empirical global tidal model developed at Ohio State University using an extended orthogonal representation of the response methods [41].
• OSUNA stands for the Oregon State University model for North of Australia, including the GBR. Developed through a representer approach [36], this model is the only regional model in this study. The assessments against coastal tide gauges have shown that this model outperforms the previous TPXO7.2 model in this region [24].
• GOT4.10 is the updated GOT99 model implemented based on an empirical harmonic analysis of satellite altimetry relative to an adopted model developed by the Goddard Space Flight Center [42]. Unlike its predecessors, e.g., 4.9, 4.8 and 4.7, in GOT 4.10, Jason-1 and Jason-2 have been used without considering T/P data.
• TPXO9 is the most recent version of the TPXO series of global ocean tide models. This model fits the Laplace tidal equations and altimetry data and has been upgraded using updated bathymetry and by assimilating more altimetry data than in previous versions [36,43].

Shallow water regions are challenging for ocean tide models [20]. Nevertheless, improvements in observations, bathymetric models and model algorithms have led to improvements in recent tidal models, not only over deep oceans but also in coastal zones [21].

3. Tidal Constituent Evaluation

The harmonic constants of nine tidal constituents (i.e., K₁, O₁, P₁, Q₁, M₂, S₂, N₂, K₂ and M₄) were computed from the SSH observations of T/P, Jason-1, Jason-2 and sea level from coastal tide gauges. These constants were compared to those of the models, revealing the strengths of the models, to assess their tidal constant estimation over the GBR.

3.1. Tidal Constants from Altimetry and Tide Gauge Data

The harmonic approach decomposes a time series of sea level observations into the contributing waves corresponding to tidal constituents. This method is used to extract major and shallow water tidal constituents from tide gauge data due to the high temporal resolution of the tide gauge data of 1 h or 0.5 h. In this study, the MATLAB Unified Tidal Analysis and Prediction (UTide) software package by Codiga [44] was employed to extract tidal constituents from tide gauge time series. Due to different instrumental and preanalysis procedures, there might be gaps in the tide gauge time series. While other available software packages require special treatments for irregularly spaced observations, UTide has the ability to analyze time series containing time gaps.

For altimetry data, another method to derive harmonic constants was used because the altimeter background noise of T/P series is within a few centimeters [45] and thus only extracts constituents with a global root mean square (RMS) above the noise level (e.g., M₂, S₂, N₂, K₁ and O₁). Therefore, the response approach was used that allows the signal-to-noise ratio problem to largely be avoided [46]. The semidiurnal (M₂), diurnal (K₁) and shallow water (M₄) constituents were used to produce the co-range and co-phase maps shown in Figure 3.

It can be seen that M₂ varies smoothly over the deepwater and shelf zones. However, in the coastal zones, large amplitudes (up to ~160 cm, Figure 3a) and rapid phase changes (0° to 360°, Figure 3b) are observed around Broad Sound. Dense offshore reefs inhibit tidal flow in the southern GBR, leading to magnification of the semidiurnal tidal constants including M₂ [6], as shown in Figure 3a,b, reaching up to 160 cm in amplitude, causing the region to be dominated by semi-diurnal tides.

The K₁ amplitude reaches 60 cm near the coast (Figure 3c) but this diurnal constituent does not undergo noticeable variations in amplitude and phase, and even in the southern GBR, the diurnal constituent (K₁) phase remains largely unaffected by dense islands. Several amphidromic points can be seen in the M₄ co-phase map (Figure 3f) and are mainly located in the coastal and shelf zones, with amplitudes reaching 5 cm.

In Figure 3, the constituents show an intense variation between latitudes −19° and −23° and between longitudes 148° and 153°. This zone includes Broad Sound, which has been the subject of previous studies due to its complex tidal regime [6]. It is more likely that the models show larger inaccuracies in the estimation of tidal constants and therefore in tidal height predictions where there is a larger tidal range [23]. In the mentioned geographic range, an intense variation in amplitude can be seen (Figure 3a,c,e), and this area is identified hereafter as the challenging zone.

3.2. Comparison of Model Tidal Constituents Against Tide Gauges and Altimetry

Nine tidal constituents of the models mentioned were bilinearly interpolated at tide gauge locations and altimetry track points. The interpolation values are hereinafter called model values. Model values were compared with values obtained from the tide gauge and altimetry observations from the combined TOPEX/Jason time series (hereinafter called validation values). The differences between the two quantities were computed for a given constituent j using the RMS error [22]:

$$RM{S}_j = {\left(\frac{1}{2L}{\displaystyle \sum _{i=1}^L{\left(H_i^m - H_i^v\right)}^2}\right)}^{1/2}$$

(2)

where L is the number of points in a zone, the superscripts m and v are the model and validation tidal constants, respectively, and H is the complex expression of tide amplitude and phase at the location as:

$$H=A(\cos (G)+\mathrm{isin}(G))$$

(3)

where A is the tidal amplitude and G the Greenwich phase of the constituent for the ith point. Considering the real part of Equation (2) gives:

$$RM{S}_j={\left(\frac{1}{2L}{\displaystyle \sum _{i=1}^L\left[{\left(A_i^m\cos \left(G_{{}_i}^m\right)-A_i^v\cos \left(G_{{}_i}^v\right)\right)}^2+{\left(A_i^m\sin \left(G_{{}_i}^m\right)-A_i^v\sin \left(G_{{}_i}^v\right)\right)}^2\right]}\right)}^{1/2}$$

(4)

In addition, to assess the performance of each model, the root sum square (RSS) was calculated using:

$$RSS={\left({\displaystyle \sum _{j=1}^{T}RM{S}_j{}^2}\right)}^{1/2}$$

(5)

where T is the total number of constituents.

RMS and RSS values were computed over the coastline, coastal, shelf and ocean zones with the smallest RMS values indicating the best model performance.

Figure 4 shows the RMS of eight major constituents averaged over the four different zones. Along the coastline, M₂, S₂, K₁, O₁ and N₂ are the most inaccurate constants estimated by the models, while, in the other zones, the most inaccurate constants estimated by the models are limited to M₂, S₂, K₁ and N₂. Three semidiurnal constituents (M₂, S₂ and N₂) are most inaccurate in all zones, which is consistent with the fact that the tidal regime of the area is mostly semidiurnal and partly mixed patterns. The OSUNA and TPXO8 models were found to be the best estimators of M₂ on the coastline with RMS values of ~13.5 cm and 13.6 cm, respectively. The HAMTIDE and EOT11a models, however, have the largest RMS values of 40 cm and 42.4 cm, respectively. FES2012 shows a better accuracy estimation for S₂ with an RMS value of 3.45 cm; OSUNA provides a better accuracy estimation for K₁ with an RMS value of 6.48 cm; and EOT11a and HAMTIDE provide better accuracy estimations for O₁ with RMS errors of 1.5 cm and 0.7 cm, respectively. Similar performance was observed for all models in estimating the P₁, Q₁, K₂ and N₂ constituents.

FES2014 shows an RMS error of ~16.8 cm for M₂, which is an improvement with respect to FES2012 (~18 cm), while for S₂ the opposite was found: with an RMS error of 3.2 for FES2012 and an RMS error of 11.2 for FES2014. The only regional model in this study, OSUNA, stands as the better estimator of four major tidal constants (M₂, S₂, K₁ and O₁ over the coastal and shelf zones (Figure 3b,c). Unfortunately, this model lacks the other four major constituents, N₂, K₂, Q₁ and P₁; these constituents are generally most accurately mapped by FES2014 over the coastal and shelf zones.

All models except for HAMTIDE include the major shallow water constituent, M₄. Several tide gauges show a presence of the principal ter-diurnal constituent, MK₃, but this has not been included in any of the models.

Table 3. The mean RMS values of differences between the model and validation values for M₄ over different zones (cm). The smallest mean RMS is indicated by bold numbers.

Zone	TPXO8	EOT11a	FES2012	OSUNA	OSU12	FES2014	GOT4.10	DTU10	TPXO9
CL	2.07	1.58	0.18	1.55	1.39	1.42	1.68	1.74	0.92
CZ	1.84	0.70	0.53	0.41	0.20	0.42	0.42	0.44	0.32
SZ	1.8	0.69	0.40	0.27	0.45	0.34	0.55	0.54	0.39
OZ	0.14	0.53	0.07	0.15	0.12	0.17	0.15	0.15	0.1

shows that the overall best estimation for M₄ occurs with models FES and OSUNA. This is expected as both models were created by the assimilation approach, which have accounted for the hydrodynamic features of the area including the bathymetry and bottom friction. However, TPXO8 falls behind the FES models in the GBR region. The purely empirical model, OSU12, shows an acceptable accuracy in calculating M₄ over all zones, particularly in the coastal zone. Except for the FES series, the other models show a similar ability in the estimation of M₄ over all four zones.

As OSUNA only has four constituents, a comparison with only the four major constituents is illustrated in

Table 4. The mean RSS values (in cm) for the four major tidal constants for the tidal models over coastline (CL), coastal zone (CZ), shelf zone (SZ) and deep ocean zone (OZ).

Zone	Coastline	Coastal Zone	Shelf Zone	Deep Ocean Zone
TPXO8	19.72	7.96	5.15	2.83
TPXO9	21.96	7.25	4.93	3.23
OSUNA	21.13	5.28	4.05	2.90
EOT11a	45.46	7.78	6.02	3.25
DTU10	23.31	6.67	5.06	3.37
GOT 4.10	30.14	8.87	6.74	3.28
FES2012	20.71	9.28	6.62	3.43
FES2014	24.43	6.60	4.62	2.67
HAMTIDE	44.78	9.69	6.05	3.14
OSU12	33.36	7.31	5.13	3.66

. This Table shows that, overall, OSUNA outperforms the other models for estimating the four major tidal constituents in the coastal and shelf zones with RSS values of ~5.2 cm and 4 cm, respectively. In the coastline and deep ocean zones, however, its RSS differs by ~1.5 cm and ~0.2 cm, respectively, from more accurate models in these corresponding zones.

To test model performance for estimating the eight major constituents in GBR, the RSS value for each model was estimated and is listed in

Table 5. The mean RSS values (in cm) for all 8 major tidal constants for the tidal models over coastline (CL), coastal zone (CZ), shelf zone (SZ) and deep ocean zone (OZ).

Zone	Coastline	Coastal Zone	Shelf Zone	Deep Ocean Zone
TPXO8	21.04	8.40	5.46	2.98
TPXO9	23.38	7.72	5.31	3.44
EOT11a	46.32	8.21	6.33	3.47
DTU10	24.52	7.13	5.37	3.61
GOT 4.10	43.14	9.34	7.12	3.49
FES2012	21.39	9.60	6.86	3.71
FES2014	23.28	7.21	5.03	2.88
HAMTIDE	45.38	10.36	6.37	3.35
OSU12	34.28	8.17	5.64	4.18

. Over the coastal, shelf and deep ocean zones, FES2014 has higher performance than its predecessor, FES2012, with an RSS below 1–2 cm. However, on the coastline it shows an RSS larger than FES2012 of ~2 cm. In the deep ocean zone, all models performed well, with similar RSS values (less than 4.5 cm, Table 5). Table 5 demonstrates that on the coastline, TPXO8 shows a superior performance with an RSS of 21.04 cm. Models FES2012, FES2014, TPXO9, OSUNA and DTU10 show similar RSS levels to TPXO8. The OSU12, GOT 4.10, HAMTIDE and EOT11a have RSS values above ~35 cm in estimating the four major tidal constants on the coastline.

The fact that the TPXO series, OSUNA and FES2014 have integrated coastal tide gauges into the modeling is possibly the reason they perform more accurately on the coastline. Although FES2012 does not include tide gauges, it uses recent bathymetric data for the coastal regions. As such, the model uses a more accurate shoreline, perhaps contributing to its better performance on the coastline. At the coastal zone, the models show similar behavior, and their RSS values differ by a maximum of ~4 cm. In the shelf zone, the RSS values are even closer, with differences among the models of just ~2 cm.

4. Model Performances in Tidal Height Prediction

Sentinel-3A SLAt time series were used to evaluate the ability of the tidal model to predict tidal heights. SLAt is the SLA that is not yet corrected for tides. Sentinel-3A is particularly interesting to use for this evaluation, as it is an independent dataset, and Sentinel-3A is employed with a new-generation synthetic aperture radar altimetry, which enables more accurate SLA data closer to the coast compared with traditional altimeters, such as used in the Jason satellites. This assessment reveals how accurate tidal constants from models (Table 5) and other features such as spatial resolution and bathymetric data (in the case of the assimilation models) contribute to tidal modeling in this region. To test the tidal prediction ability of models, the following procedure was considered.

(1) SLAt values computed at the tide gauge stations (Equation (1)) were used for comparison at the coastline. In this section, one year of hourly sea level time series from tide gauges were used.

(2) For each Sentinel-3A along-track location, a time series of SLAt was detided by the time series of tidal heights (TH) from the individual tide models. The associated software with models was used to predict tidal heights. In the coastal, shelf and deep ocean zones, 40 cycles of the Sentinel-3A mission were used. The resulting detided time series of SLA were defined as:

$$SLA=SLA{t}^{}-TH$$

(6)

The RMS value was computed for the estimated SLA at each Sentinel-3A along-track position (Equation (7)).

$$RM{S}_{}={\left[\frac{1}{N}{\displaystyle \sum _{n=1}^{N}SL{A}_n^2}\right]}^{1/2}$$

(7)

where N is the number of observations.

The SLA includes tidal height residuals, unmodeled atmospheric effects that have remained due to errors in atmospheric corrections, surface currents and other oceanographic signals. In addition, the tidal residuals of the other signals were assumed to be the same for each time series.

Additionally, in this case, the RMS values were averaged over the single zones (

Table 6. Mean RMS values (in cm) of the SLA, calculated by detiding SLA using different ocean tidal models, based on one year of hourly tide gauge observations in coastline and three years of Sentinel-3A data (40 cycles) over other zones, in the GBR. The smallest mean RMS values are indicated by bold numbers.

Model	Coastline	Coastal Zone	Shelf Zone	Deep Ocean Zone	RSS
TPXO8	31.85	13.82	9.11	7.90	36.73
TPXO9	33.26	15.68	9.99	7.96	38.92
OSUNA	32.82	20.72	17.66	17.9	45.19
EOT11a	36.06	23.38	16.18	10.61	47.10
DTU10	25.78	13.30	9.48	7.82	31.47
GOT 4.10	36.33	14.52	11.48	7.73	41.49
FES2012	26.62	10.94	7.71	7.64	30.81
FES2014	26.05	16.33	12.99	7.46	34.16
HAMTIDE	41.89	17.19	13.13	8.88	47.96
OSU12	37.35	25.47	17.95	12.64	50.17

). The models with more accurate tidal height predictions show lower RMS values over Sentinel-3A along-track and tide gauge locations. Table 6 shows the statistics for the mean RMS values of different models’ SLA. As expected, the models show lower accuracies over shallow regions of the coastline and the coastal zone, reflected in higher RMS values.

In general, at the coastline, RMS values range between 26 and ~41.9 cm. HAMTIDE shows the maximum RMS of ~41.9 cm, while the DTU10 and FES models have the minimum RMS values (~26–27 cm). In the coastal and shelf zones, FES2012 compares more favorably than others with RMS values of ~10.9 cm and ~7.7 cm, respectively, while, in the deep ocean zone, FES2014 is more favorable with an RMS of ~7.4 cm. In the coastal and shelf zones OSU12, and in the deep ocean zone OSUNA stand behind the other models with RMS values of ~25.4 cm, ~17.95 cm and ~17.9 cm, respectively. Overall, FES2012 with an RSS value of ~30.81 cm outperforms the other models followed by DTU10 with an RSS value of ~31.5. The pure empirical model, OSU12, has a lower precision than others with an RSS value of ~50.2 cm.

Figure 5 shows the variations in RMS values over the GBR. In general, FES2012 outperforms the other models in this region. As expected, the regional model, OSUNA, is the most inaccurate model in terms of tidal height prediction ability over the shelf and deep ocean zones because of the few constituents available (Table 4). In fact, FES models, which include more than 30 major and shallow water tidal constants, provide more accurate tidal height predictions. According to Figure 5h, FES2012 is the best model for predicting tidal heights. Although FES2014 includes more constituents than FES2012, the predicted tidal heights are less accurate below −20° latitude (Figure 5i), where the RMS values exceed 40 cm. In the challenging zone, FES2014 shows a noticeably large average RMS (~21.2 cm), while its predecessor, FES2012, in contrast, performs with average RMS values of ~8.2 cm with large levels only appearing along the coastline.

As shown in Figure 5a,b,f, the DTU10, TPXO9 and TPXO8 models show large RMS values over the challenging zone. Similar behavior can be observed for the other models. In the case of GOT4.10, this may be due to its lower spatial resolution. The HAMTIDE model is more accurate than EOT11a and this may be because HAMTIDE is based on EOT11a but also takes into account hydrodynamic features that may be helpful for characterizing this area. The EOT11a model mainly shows inaccuracies close to the coastline, and very large RMS values appear in the challenging zone. In Figure 5a,b, TPXO9 shows larger RMS values over the challenging zone than its predecessor, TPXO8. Over the other zones however, both models show similar performance.

Figure 5, shows that large RMS uncertainties appear near the coast as well as in Broad Sound. The profile of SLA RMS values was also plotted with respect to the bottom topography and the distance to the coast to investigate their variations according to the tidal model used (Figure 6). For this experiment, pass 731 of the Sentinel-3A mission was used (cf. red pass in Figure 2). Starting from the area around Broad Sound, pass 731 intersects the challenging zone below −20° latitude and crosses over the coral reefs towards the northeast. This passes over highly varying bottom topography with bathymetry varying from 1 m to 50 m before reaching the continental shelf at a distance of ~120 km from the coastline. As such, this pass can be used to investigate how the models behave in this complex zone.

Major differences between models appear within ~100 km of the coastline, which has sudden bottom topography variations due to the reefs. In detail, all RMS profiles fluctuate with a maximum of ~94.1 cm for EOT11a and a minimum of 5.21 cm for FES2012 within ~3 km of the coast, showing strong sensitivity to the fluctuating bottom topography. The EOT11a model shows problematic performance with a mean RMS value of 73.2 cm, followed by FES2014 (41.5 cm) within 15 km of the coast (coastal zone) along Sentinel-3A pass 731. The FES2012 model shows the best performance in compliance with Figure 5h and Table 6. Other models have mean RMS values of less than 20 cm. This RMS pattern remains unchanged for a distance of ~30 km to 80 km from the coastline, although the RMS values show relatively smooth variations with respect to less variable bathymetry. As expected, all RMS values decrease over the open ocean (beyond 120 km from the coast). It was observed that five models (FES2012, DTU10, GOT4.10, TPXO9 and TPXO8) have much smaller RMS values than the other models along pass 731, especially within a distance of 80 km from the coastline. The intense variations in the RMS profiles over the highly fluctuating bathymetry (Figure 6) confirm the significance of bottom topography variations and their contribution to model accuracy regarding tidal constituent estimation and prediction. Since all models (empirical or assimilation) used altimetry observations and there are many islands within the study area, the large RMS values of the SLA may partly be affected by altimetry errors from retracking distorted waveforms close to the islands and the coastal zone.

The RMS profiles over distances of 80 km to 120 km vary smoothly, although there are also sudden depth changes in bottom topography. As such, a combination of factors such as proximity to the coastline or to islands and sudden changes in bottom topography contributes to inaccuracies in the tidal model prediction ability over this area.

5. Discussion

The performance of recent tidal models in tidal estimation and prediction over the GBR region was evaluated. Overall, the model with better performance was found to be FES2012. This model shows height predictions over the GBR region to an accuracy of 10.9 cm in the coastal zone and provides a reasonable spatial resolution that captures small tidal waves, an accurate bathymetry model that effectively illustrates the sudden variations in depth over the challenging zone and includes many influential major and shallow water constituents (Table 5 and Table 6, and Figure 5e and Figure 6).

The latest version of the FES models, FES2014, estimates constituents at along-track points with higher accuracy than does FES2012 over the coastal, shelf and deep ocean zones, with RSS differences of ~2 cm. However, this model shows larger RMS values in the challenging zone compared to FES2012. According to Cancet et al [47], this is a challenging zone for FES2014 that also shows problems in predicting tidal currents.

DTU10 and EOT11a share the same base model (FES2004) in a remove and restore procedure. However, different approaches such as the response method and harmonic analysis are used to estimate corrections to tidal constants. This may be the cause of ~1 cm of RSS difference in tidal constant estimation on T/P along-track positions located in the coastal and shelf zones (Table 5). DTU10 uses a depth-dependent method to spread tidal constants onto a regular grid which in turn can play a role in providing more accurate results over the highly variable depths of the challenging zone. Therefore, by accounting for the impact of sudden bathymetry variations on tidal constants, DTU10 compares more favorably than does EOT11a over the GBR and the challenging zone (Figure 5c,e).

The regional model, OSUNA, accurately estimates four major tidal constants (Figure 4). However, the low number of major constituents used in this model is not sufficient to accurately predict tidal height. This flaw is shown by the numerical values in Table 6 and in Figure 5b, in which the highly variable RMS values are geographically depicted. TPXO9 shows lower performance than its precursor, TPXO8, despite its finer bathymetry. Although using more accurate bathymetry data in the assimilation models is a precursor to develop more efficient models, it is not the only factor. In the OSU models, representer locations, bottom friction and other datasets such as current meters contribute to efficiency of the model. In Figure 5 and Table 5 and Table 6, it can be seen that the addition of new datasets and improvements to TPXO8 has not resulted in higher performance for TPXO9 in this area. Over coastal, shelf and ocean zones, HAMTIDE has been constrained by the tidal constants estimated from EOT11a [40]. As such, in Table 6, it can be seen that accuracy of the empirical model EOT11a has been improved.

It is noteworthy that the RMS values of SLA used to evaluate the tidal height predictions of the tidal models includes tidal residuals, unmodeled surface circulations and high frequency atmospheric effects. The surface circulation patterns in the region are the results of interactions among surface currents and highly variable bathymetry, coral reefs and the dense presence of islands, particularly in the southern GBR. As such, the sea level in the area of study is also subject to variations due to large scale currents such as the East Australian Current (EAC) and the North Queensland Current (NQC), as well as small scale currents. Although the SSHs have been corrected for DAC, the minor and high frequency atmospheric effects undoubtedly generate sea level variations, with greater effects closer to the coast. Consequently, the large RMS values shown in Table 6 could be due to a combination of these possible factors.

Per Figure 6, models show intense RMS variations within an 80 km distance from the coast where the bathymetry has sudden fluctuations. In addition, according to Figure 5, the larger RMS values are mainly confined to the challenging zone. The challenging zone, with a maximum depth of 50 m, represents the most variable area in terms of bathymetry due to the existence of small islands. Given that all assimilation and empirical models considered in this study have used altimetry SSH which is more likely to be affected in the vicinity of dense islands, the models generally show greater inefficiency over the challenging zone in comparison to other areas. In proximity to the coastline and the islands in the southern part of the GBR, where altimetry observations are suspected to be affected by distorted waveforms, the sudden variations in bathymetry are more destructive to model predictive abilities compared to the areas that are more distant from the coastline and islands (Figure 6).

It was found that a purely empirical model (OSU12) that neglects existing complexities is not able to capture the real image of tides in this area. As such, an assimilation model that accounts for highly variable bathymetry and bottom friction is the best approach to model the tides in this region. Drag coefficient variations, of interest over sand bed oceans, are subject to large fluctuations over areas with different seafloor textures, such as coral reefs. This parameter may be not yet well described in existing models, and it would be worthwhile to further investigate together with assimilation of more accurate bathymetry for the GBR region, a study of which is currently ongoing. In addition, the retracked altimetry data within 10 km from the coastline contribute to more efficient tidal models that are able to reflect the complexities of this phenomenon over the GBR.