Defining humpback whale (Megaptera novaeangliae) potential distribution in the Great Barrier Reef Marine Park: a two-way approach

Study site

The GBRMP is located on the east coast of Australia (Fig. 1) and contains one of the world’s richest and most complex ecosystems. It has an average depth of 35 m in its coastal waters, whereas on outer reefs, continental slopes extend to depths of more than 2000 m (Great Barrier Reef Marine Park Authority 2019). The GBRMP has a tropical climate, with average temperatures in summer between 24°C and 33°C and in winter between 14°C and 26°C. The region has two distinct seasons, a dry, warm season and a hot, wet season. The Great Barrier Reef covers approximately 388,000 km², with 99% of the area designated as GBRMP (Great Barrier Reef Marine Park Authority 2019). From June through November, the region is an important HW breeding area (Simmons and Marsh 1986; Smith et al. 2012).

Fig. 1.

Location of the study area (GBRMP) within Australian waters, showing sightings of HWs (yellow circles) used for modelling. The inset shows Australia with the location of the GBRMP displayed in yellow.

Whale records

HW sightings were retrieved from Eye on the Reef (EOR), Atlas of Living Australia (ALA), Ocean Biodiversity Information System (OBIS), and the National Marine Mammal Database (NMMD). The largest number of records obtained for the study were from EOR, consisting of 3640 observations collected between 2000 and 2020. EOR is a monitoring and assessment program in the GBRMP, where staff from tourist dive boats collected standardised biological information regarding the reef’s health at frequently visited sites. This includes recording opportunistic HW sightings during each trip. ALA, OBIS and NMMD are open-access collaborative databases on biodiversity. These databases employ a combination of methods to collect data from various sources, including citizen science programs, ship-based and aerial surveys, partnerships with organisations, and data digitisation. The associated data from these (88, 189, and 19 sightings, respectively) are presence-only opportunistic records collected non-systematically from 2000 to 2020. Although only limited information is available regarding their research efforts on HWs, collected data have undergone quality checks to ensure accuracy and reliability. Quality checks include location, time and date validation, data cleaning and filtering, cross-checking with other databases, and expert review to ensure data accuracy.

The datasets utilised in this study, in contrast to the approach by Smith et al. (2012), extend beyond the peak breeding months of July and August. This broader temporal scope aligns with research indicating potential changes in the timing of both migration and breeding activities (Franklin et al. 2011; Ramp et al. 2015; Meynecke et al. 2017; Torre-Williams et al. 2019; Gosby et al. 2022).

Environmental variables

Salinity, SST and chlorophyll a values were obtained from BioOracle (Tyberghein et al. 2012). The long-term averages for these data span a time period of 2000–2014 at a spatial resolution of approximately 8.3 km. Bathymetry was obtained from eAtlas (Beaman 2020) at a raster resolution of 100 m. Pixels with positive bathymetry values, corresponding to islands between the GBRMP and part of the coastline, were removed from the raster.

Distance to the coast and distance to the reef were derived using ArcGIS Pro ver. 2.9.0 software geoprocessing tools. For distance to the coast, each pixel in the raster represents the distance in metres from the centre of that pixel to the coastline of the GBRMP (excluding any islands). Geoprocessing tools were utilised to calculate the distance of each point within the study domain to the coast and convert these distance values into a raster format. Similarly, distance to the nearest reef was computed using the same process, measuring the distance from each point to the closest reef rather than the GBRMP coast.

In this study, ‘Coastal Waters’ or ‘inshore’ areas adhere to the definition provided by Geoscience Australia (2023), based on the United Nations Convention on the Law of the Sea (UNCLOS). This approach provides a consistent framework for comparison between datasets. Under this definition, ‘Coastal Waters’ areas extend from the territorial sea baseline up to a line 3 nautical miles seaward.

All grids were resampled to a consistent 0.085° (~8.5 km) grid resolution and reprojected to the standard WGS84, EPSG 4326 geographic projection in ArcMap, ensuring uniform geographic bounds and cell size for all variables.

Bias correction

Two different bias-correcting methods were applied to account for spatial clustering of HW sightings data (Fig. 1). First, data were spatially rarefied at 10-km resolution using the ‘spThin’ R package (Aiello-Lammens et al. 2015; Duengen et al. 2022). This resulted in a total of 2584 (92.9%) records being removed from the 2000–2014 dataset, leaving 197 records (Model C; Table 1). Second, we employed the ArcGIS ‘Sample by buffered MCP’ tool for background selection. This approach involved selecting background points within a maximum radial distance of 5 km from known occurrences, with a cap of 10,000 points. This method aligns with established practices in the field (Phillips et al. 2009; Barbet-Massin et al. 2012; Brown 2014).

Maxent tuning

Maxent’s default parameters were used for Models A–H (Table 1). Subsequently, the best-fit models from each time span were optimised through tuning to enhance model performance. The best-fit Maxent models were selected based on the Bayesian Information Criterion (BIC) calculated with ENMTools (Warren et al. 2010). BIC was chosen for model selection because it penalises model complexity more strongly than other criteria (Burnham and Anderson 2002; Neath and Cavanaugh 2012; Kumar et al. 2015). The ENMeval package in R (Muscarella et al. 2014) was used to conduct the tuning process. For the best-fit models from each time span, 24 combinations were evaluated: six feature classes (L/LQ/H/LQH/LQHP/LQHPT; where ‘L’ linear, ‘Q’ quadratic, ‘H’ hinge, ‘P’ product transformation and ‘T’ threshold), and four Regularisation Multiplier (RM) values (from 0.5 to 3.5, with one unit increment). Using the two combinations that achieved the highest Akaike’s information criterion (AIC), the tuning process resulted in creating two Maxent models, Models I and J.

Model evaluation

Response curves were used to show how each variable affects the Maxent models’ overall predictions. They demonstrate how prediction probabilities of presence vary as each environmental variable is modified, using only one variable at a time. The jackknife test was used to assess the contribution of each variable in the models, providing insights into their explanatory power and unique information for the species’ potential distribution.

Maxent models were assessed using threshold-independent and threshold-dependent metrics to evaluate their performance. The area under the curve (AUC) of the receiver operating characteristic (ROC) curve served as a threshold-independent measure, quantifying the model’s discriminatory power. AUC values range from 0 to 1, where 0.5 signifies random prediction, and 1 indicates perfect discrimination between suitable and non-suitable habitats. Two metrics were employed for threshold-dependent evaluation: omission rate (OR) at the minimum training presence threshold and OR at the 10% training presence threshold. The expected value of the test OR at the minimum training presence threshold is 0, whereas at the 10% training presence threshold, it is 0.10. Higher-than-expected OR values indicate inadequate model performance (Boria et al. 2014). The Maxent algorithm calculates the OR in training and test data, representing the percentage of test locations predicted as unsuitable for the species (Phillips 2005).

To validate the representativeness of our dataset and comprehensive environmental coverage across the GBRMP compared to Smith et al. (2012), we conducted a niche overlap analysis (Lysy et al. 2014). The analysis was performed using the ‘nicheROVER’ R package. The variables for the analyses were distance to the shore, distance to the reef, SST, bathymetry, and salinity. The datasets used were the BPC sightings employed by Smith et al. (2012) compared to the dataset used for our best-performing model.

Model development

Data processing and statistical analysis were conducted using the following software versions: QGIS 3.22.5 (QGIS Development Team 2022), ArcGIS Pro 2.9.0 (Esri 2021), and R 4.2.0 (R Core Team 2022). Predictive models of HW distribution in the GBRMP were developed using the software Maxent 3.4.4 (Phillips et al. 2021). The first models (Models A and B) were developed using the full HW sightings dataset from 2000–2014 and 2000–2020 and used as the baseline models with no post-processing data (Table 1). Models C to J, our post-processed models, were developed to improve the performance of Models A and B, which would result in the creation of our best-fit model (Fig. 2). A comparison was made between two sets of models (Models A–B; C–D; E–F; G–H; I–J) to ensure the reliability of the 2000–2020 models despite the environmental predictors not aligning with the longer period of whale sightings. This approach was necessary due to the limited availability of environmental data in BioOracle, which only covered the years 2000–2014. The selection of the 2000–2014 period for the environmental data by BioOracle was based on the availability of satellite data used to generate the environmental models, ensuring the use of the most reliable and up-to-date information for that specific timeframe. Multi-collinearity among environmental variables was tested by calculating the Kendall correlation rank coefficient (Table 2). Correlated variables were identified using a threshold value of ±0.5 (Dormann et al. 2013). The variable that contributed the least information to the SDM according to the jackknife test was excluded for subsequent models’ development.

Fig. 2.

Methodological flow chart showing the sequential steps followed in this study: Maxent models were initially developed using 2000–2014 data, later expanded with 2015–2020 sightings. Model C, refined through bias corrections and tuning, was employed to reassess HW sighting data from Smith et al. (2012), validating Model C’s accuracy and facilitating a final comparison with Smith et al.’s model.

The best-fit model was compared with the model prediction of environmental suitability developed by Smith et al. (2012). Differences in response curves and variable contributions were analysed. Additionally, our model results showcased the most suitable areas for HWs within the GBRMP, highlighting these habitats through maps. This approach enabled the evaluation of habitat suitability congruence between the two models for HWs. To demonstrate the efficacy of our approach, Model X (Table 1) was developed by utilising an approximation of the BPC dataset (90 HW sightings collected between 2003 and 2007) and incorporating the environmental predictors and parameters from our best-fitting model. Locations of the HW sightings in the dataset used in Smith et al.’s (2012) research were estimated using georeferencing in QGIS before using them as input data in Maxent.

Whale records

Observations followed the expected seasonal trend of HW presence in the GBRMP, with most observations (80%) occurring in July and August (Fig. 3). The monthly distribution of HW sightings post-filtering within the GBRMP (Supplementary Fig. S1) also depicts the peak of observations during July and August.

Fig. 3.

Monthly distribution of HW sightings within the GBRMP, comparing the periods 2000–2014 (blue) and 2000–2020 (orange).

Model development

According to the Kendall correlation rank test, salinity, SST, and chlorophyll a were identified as correlated variables (Table 2). Based on the jackknife test, chlorophyll a contributed the least information to the full model and was excluded from subsequent models.

The niche overlap analysis revealed a significant 96% congruence (α = 0.99) between the environmental niches defined by our dataset and those identified in Smith et al. (2012).

Both Models A and B suggested highly suitable areas (>0.7) were primarily concentrated near the coast, between 19°S and 21°S (Fig. 4). The average test AUC (Table 3) for both models exceeded 0.8, indicating the models are reliable at discriminating between suitable and non-suitable habitats. SST and distance to the coast contributed most to the overall performance of Models A and B, collectively explaining approximately 80% of the variation in species distribution (Table 4). The variables that contributed the least were bathymetry for Model A and distance to the reef for Model B. For both models, the variable with the most predictive power (highest AUC) when used in isolation was distance to the coast; therefore, it was the most significant driver for each model (Fig. S2). The variable with the least predictive power was distance to the reef for both models. The variable that decreases the AUC the most when omitted is SST for Model A and distance to the coast for Model B; therefore, they appear to have the most information not represented in the other variables. The response curves for both models A and B are similar (Fig. S3). Distance to the coast and SST values that correspond to high environmental suitability (>7) changed by no more than 0.5 units on average.

Fig. 4.

Environmental suitability for HWs in the GBRMP between (a) 2000–2014 (Model A) and (b) 2000–2020 (Model B). The colour scale denotes habitat suitability, ranging from 0 (light blue) to 1 (dark red).

In total, 10 Maxent models, with varying feature classes and RM, were evaluated for predicting HW distribution using different combinations of bias-correction methods. All models performed better than random, with training and test AUC values exceeding 0.7 (Table 3). Moreover, these models demonstrated low OR at the minimum predicted threshold, with values of 0 in the training data and below 0.02 in the test data. Additionally, at the 10th percentile training presence threshold, the models showcased OR values below 0.09 in the training data and below 0.16 in the test data. Model F had the highest BIC value and lowest rank. After bias correction, each method (spatial rarefying, implementation of a bias file, and both methods combined) resulted in a slightly lower testing AUC and a broader range of suitable habitats compared to models with no bias correction (Table 3). Evaluations of RM and feature class combinations for Maxent tuning are presented in Tables S1 and S2.

The analysis of relative contributions revealed that SST and distance to the coast emerged as the two key variables influencing the performance of Models C to J. These variables, on average, accounted for approximately 75% of the variation in species distribution across these models (Table 5). The jackknife plot of variable importance indicated that distance to the coast, when used alone, was the environmental predictor that had the greatest influence on Models C to J (Fig. S4). The variable that decreased the models’ AUC the most when excluded was also distance to the coast. On the other hand, distance to the reef had the lowest impact on models D, J, and I. In models C, E, G and H, the lowest impact was for salinity, and for Model F, it was SST.

Model results for the two different periods, 2000–2014 and 2000–2020, showed similar results (Fig. S5).

Comparison between the best-performing Maxent model and the existing model from Smith et al. (2012)

Model C emerged as the best-fit model (BIC value: 3081.2). It included the feature class LQH with an RM value of 0.5 (Table 1). Model C suggested the highest values of environmental suitability (0.7–0.99) between 15°S and 21°S (Fig. 5). In contrast, Smith et al. (2012) identified two core areas as highly suitable for HW based on their Maxent model. The first area is located approximately 100 km east of Proserpine, south of Mackay, and the second area is offshore of the Capricorn and Bunker groups of islands and reefs (Fig. 5). The response curves for Model C showed a preference of HW for waters between 7 and 33 m deep, with the highest probability at −11 m. In contrast, Smith et al.’s (2012) identified the highest values of environmental suitability (>0.7) at depths between 30 and 58 m, with the highest likelihood at −49 m (Fig. 6). Similarly, the marginal response curve for sea surface temperature (SST) differed between the models. Model C exhibited a range of 25.7–26.4°C with the highest value of 26.2°C, whereas Smith et al. (2012) reported a range of 21–23°C with the highest probability at 21.8°C. Additionally, Model C suggested HWs showed a preference for areas between 0.24 and 37 km from the coast, whereas Smith et al.’s (2012) identified a preference for areas between 60 and 110 km from the coast.

Fig. 5.

(a) Model prediction of average environmental suitability of the best-fitting model, Model C. (b) Model prediction of average environmental suitability between 2003 and 2007 developed in research by Smith et al. (2012). Colour scale denotes habitat suitability ranging from 0 (light blue) to 1 (dark red). GBR, Great Barrier Reef.

Fig. 6.

Response curves depicting the relationship between the predicted probability of HW occurrence and three environmental variables – SST, bathymetry, and distance to the coast. (a) Response curves for Model C developed in this study. (b) Response curves for Smith et al.’s (2012) model.

Results from the jackknife test showed that the variable that contributed the most to Smith et al.’s (2012) model was SST, followed by water depth (Fig. 7). In contrast, in Model C distance to the coast was the most important variable, followed by SST.

Fig. 7.

Results of the jackknife test of variable importance to AUC for (a) Model C and (b) Model by Smith et al. (2012). Blue/white bars represent the AUC values when each model uses only a single variable. Turquoise/black bars depict the AUC values when each model removes a single variable from the general model. Red/grey bars represent the average AUC value across all test iterations for each model.

As further validation, Model X was developed using the model settings and environmental variables (SST, bathymetry, distance to the coast, distance to the reef, and salinity) that yielded the best fit (Model C), along with an approximation of the BPC dataset used in the study by Smith et al. (2012). Model X gave approximately the same potential habitat suitability distribution (Fig. S6) and the relationship between the predicted probability of occurrence and the environmental variables (Fig. S7) as the model by Smith et al. (2012). In both models, bathymetry and SST were identified as the most important drivers according to the jackknife test (Fig. S8).

Species distribution modelling

We developed the best-fit Maxent model for HW distribution in the GBRMP based on citizen science data using a range of bias-correction methods. We observed that areas with the highest habitat suitability for HWs were concentrated between 15 and 21°S, predominantly in the north-eastern coastal region of the GBRMP. Environmental characteristics of suitable areas suggest that high habitat suitability for HWs generally occurs in shallow (−7 m and −33 m) and warm waters (25.7–26.4°C) that are nearshore (0.24–37 km from the coast). These results are consistent with previous research, where HW showed a preference for shallow (Derville et al. 2019; Chou et al. 2020) and warm (Rasmussen et al. 2007) habitats in their breeding grounds. HWs on breeding grounds are often in nearshore waters (Thums et al. 2018; Chou et al. 2020), as we found for HWs from subpopulation E1 in the present study.

Furthermore, our findings revealed that implementing bias-correction methods improved the complexity of the SDMs, but using the spatial rarefying approach resulted in the best-fitting model. Through this process, we generated further distribution models for HW in the GBRMP and contributed to the assessment/validation of such modelling approaches. Our data indicated clustering, which is a common challenge when dealing with opportunistic sightings. Spatial rarefying removed >90% of data in datasets from both time spans. Despite the reduction in data, spatial rarefying led to improved models and produced the best-fit model (Model C). Thus, spatial rarefying is an adequately robust correction method that helps to enhance the accuracy and reliability of SDMs based on opportunistic sightings (Kramer-Schadt et al. 2013; Boria et al. 2014; Fourcade et al. 2014). However, applying bias corrections in models C – J decreased AUC scores compared to the models with no bias correction, A and B, possibly because biased data produce an elevated assessment of model performance because of the overfitting of spatially autocorrelated data (Radosavljevic and Anderson 2014). Overfitted models often generate overly complex models, which are likely to have poor predictive ability for unsampled cells and to make predictions that fit too closely to known occurrences (Peterson et al. 2011). Consequently, the observed decrease in model performance following bias correction could indicate an improvement in addressing the biases. In addition, models with bias corrections showed a more widely spread habitat suitability distribution of HW than models with no bias correction, which is also an indicator of SDM improvement.

Comparison to previous research

We compared our model output to an existing model (Smith et al. 2012), revealing a lack of fit between the two models. By utilising the model settings and environmental variables that yielded the best fit, in conjunction with an approximation of the sighting dataset used in the study by Smith et al. (2012), we constructed a model that validated our approach. Our results suggested the significance and contribution of environmental factors to a model vary when utilising different datasets.

After establishing confidence in the Maxent settings used for Model C, we compared its results with the SDM developed by Smith et al. (2012), which utilised different opportunistic data. Our findings highlight a distinct preference for the eastern coastal region of the GBRMP in Model C, whereas Smith et al. (2012) observed the highest habitat preference near the reef itself. The disparity could be attributed to the distinct spatial coverage of whale sightings by the two datasets (Fig. S9). The limitations of whale-watching excursions, where vessels have limited search time and tend to cover the same area on each trip, are likely to have led to a spatial bias in our dataset. Vessels prioritise sightings of whales closest to the harbour based on previous observations, potentially resulting in duplicated counts and clustering of sightings around nearshore areas. In comparison, BPC aerial surveillance data predominantly captures offshore reef regions.

Furthermore, the variation in temporal coverage of the data may also contribute to the discrepancy. Model C encompassed all months from 2000 to 2014, whereas Smith et al. (2012) focused solely on observations from July and August between 2003 and 2007. This limited temporal range might have resulted in fewer recorded HW sightings along the coast, potentially underestimating the suitability of this zone as a preferred habitat. An additional factor to consider is the inshore movement and increased coastal dependence of HWs, especially among calf groups, observed later in the breeding season (Smith et al. 2021). Model C sighting data covering the full breeding season provides a different perspective on HW habitat preferences, highlighting the greater coastal dependency and demonstrating the importance of considering temporal coverage when assessing HW distribution, particularly in relation to their breeding behaviours and spatial preferences.

Still, results from Model X provided evidence for the validity of our modelling techniques by predicting a remarkably similar potential habitat suitability distribution to Smith et al.’s (2012) model when using the same data, and the response curves of both models exhibited similar relationships among variables. Importantly, the jackknife test identified bathymetry and SST as the most influential drivers of the species distribution. That means the same relationship between HWs and their environment is retained when using the BPC dataset along with the environmental variables and the Maxent parameters selected for Model C. The main difference between the two models lies in the presence of larger areas with high habitat suitability in the Maxent model, indicated by the colour red in Model X (Fig. S6). This discrepancy can be attributed to the distinct set of variables employed in each model. Model X employed environmental predictors such as SST, salinity, bathymetry, distance to the coast, and distance to the reef, which aligns closely with the variables used by Smith et al. (2012). However, whereas we incorporated salinity in our model, Smith et al. opted to include seafloor slope in their analysis. This suggests that the patterns of habitat suitability for HWs are influenced by variable selection. Although salinity did not emerge as a primary driver, its presence in the model contributes to a better understanding of the marine environment by potentially capturing slight variations in habitat preferences that are not solely dictated by temperature and depth. Our decision to include salinity, given its expected significance in coastal regions, diverges from Smith et al.’s model, which excluded this variable, likely due to most of the BPC sightings being offshore. Furthermore, Smith et al.’s findings on the homogeneity of the seafloor slope and the absence of clear patterns in the response curves suggest a minimal impact on HW distribution in the GBRMP (Smith et al. 2021). Given that slope is a secondary metric derived from bathymetry (Dolan 2012), its value is limited without significant depth variations, such as those near continental shelves. This rationale guided our preference for bathymetry use over derivatives like slope or rugosity.

The high percentage of niche overlap (95%), supported by whale sightings in these areas, suggests that nearly all the environmental conditions found in the dataset employed by Smith et al. (2012) are also present in our dataset. This finding emphasises that the environmental diversity of the GBRMP is adequately represented in Model C.

Future studies, outlook, and improvements

This study demonstrates the benefits and limitations of using opportunistic data to establish habitat use by HWs. Opportunistic data are collected by non-experts, often without a clear scientific objective. Although these data sources can provide valuable information on species occurrence and environmental conditions, inadequate or biased data can result in inaccurate model predictions. For instance, the data used herein were biased towards the coast, causing our models to overestimate or underestimate the habitat suitability in other areas. Opportunistic data may be incomplete or contain errors, which can further reduce the accuracy of predictive SDMs. For example, we encountered HW sightings falling on land (incorrect longitude/latitude), duplicate data and false entries (improper data logging), and not all errors can be removed. It is further important to recognise that Maxent modelling itself also has inherent limitations. Although capable of handling opportunistic sightings and offering options to address bias through tuning, filtering, or other methods, Maxent ideally requires systematic sampling of the entire sampling domain. In many cases, including for HW sightings data, this is not the case, thereby limiting the full predictive powers of the modelling technique. Especially in the absence of sufficient spatial coverage, incorporating scientific survey data alongside citizen science data stands out as the most effective means to obtain models that closely reflect reality. This dual approach can provide a more comprehensive understanding of habitat suitability, addressing potential biases and enhancing the overall reliability of species distribution models.

Despite their limitations and biases, opportunistic data offer a wealth of information on species occurrence. With innovative approaches and careful consideration, researchers can use citizen science data to better understand and protect HWs. This collaborative approach holds immense promise for the continued conservation of HWs and their critical habitats. The findings of this study have important implications for the conservation of HWs and their habitats by shedding light on significant areas that enhance existing knowledge and support conservation initiatives in the GBRMP. For instance, the Whitsundays (part of the GBRMP) are home to HWs in the breeding season and are poised to be declared a World Whale Heritage area due to their significance (Roulston 2023). This designation highlights the need for comprehensive studies and conservation efforts in this area.