Finding lobsters: investigating a period of unusually low settlement of Panulirus cygnus by using larval dispersal modelling

Keywords: Capes Current, fisheries, larval dispersal, Leeuwin Current, Panulirus cygnus, particle tracking, puerulus, settlement patterns, western rock lobster.

Lagrangian analysis of ocean circulation and pathways is an essential tool for understanding the dynamic motions of seawater and the nutrients, heat and biological matter within it (van Sebille et al. 2018; Hawkins et al. 2019; Jones et al. 2020). From a practical and theoretical perspective, it forms the backbone of understanding dispersal throughout the ocean at local or vast spatial and temporal scales. For example, the distribution of marine species during their early larval stages can have varying consequences through high mortality because of current movements away from recruitment locations (Chiswell and Booth 2008) or survival through retention in eddy systems (Hill et al. 1996; Cetina-Heredia et al. 2019) and connecting populations (Chiswell et al. 2003; Cowen and Sponaugle 2009; Nolasco et al. 2018; Pires et al. 2020). Furthermore, understanding these mechanisms for a single species influences their spatial management scales and their likely responses to environmental disturbances or shifts (Kvile et al. 2018; Clavel-Henry et al. 2020).

In the life cycle of spiny lobster species, the pelagic larval stages are of utmost importance to their survival (Fig. 1). The length of this offshore journey varies from several months (Palinurus elephas) to 2 years (Jasus edwardsii; Booth and Phillips 1994). Those with long pelagic larval durations can be transported up to thousands of kilometres offshore. However, behavioural changes at different development stages can control the dispersal, assisting their return to preferential juvenile habitats (Butler et al. 2011). These behaviour patterns include vertical migration regulated by light changes. During the day, larvae (phyllosoma) sink to darker depths and, by night, return to prey-rich surface waters, defined as diurnal vertical migration (DVM). These depths during day and night also can change with age, known as ontogenetic vertical migration (OVM), and vary with the moon phases (Rimmer and Phillips 1979). Phyllosomas are within offshore current systems throughout these migrations (Fig. 1). They build up their lipid reserves to undergo their final metamorphosis into a puerulus (Wilkin and Jeffs 2011; Stanley et al. 2015). At this stage, they have developed swimming abilities. In addition, they are nektonic (non-feeding; Jeffs et al. 2005; Espinosa-Magaña et al. 2018), where their biochemical energy reserves assist their transport across continental shelves and towards reef or seagrass habitats (Phillips et al. 1978).

Fig. 1.  A schematic diagram representing the three major early life stages of P. cygnus and their approximate locations off the western coast of Australia. This includes the approximate age and timing over the calendar year.

For the spiny lobster species endemic to the western coast of Australia, the western rock lobster, Panulirus cygnus (P. cygnus), the pelagic larval duration begins at hatching between September and February every year (Fig. 1). Phyllosoma are first immersed in the ocean between 40 and 100 m (Melville-Smith et al. 2009) and, subsequently, driven off the continental shelf by alongshore winds, through Ekman transport, and the region’s dominating eastern boundary current, the Leeuwin Current (LC; Phillips et al. 1981; Feng et al. 2011; Wang et al. 2015). The LC eddies retain the larvae as zooplankton while they feed on soft-bodied prey with a preference for chaetognaths (Saunders et al. 2012; Säwström et al. 2014). They undergo a series of moults with nine distinct life stages, with OVM over 8–10 months in the south-eastern Indian Ocean (Braine et al. 1979). Then, metamorphosis to puerulus occurs before crossing the continental shelf over approximately a month (Fig. 1; Chittleborough and Thomas 1969; Feng et al. 2010, 2011). The puerulus settlement season, whereby these recruits enter the fishery, is from May to April each year (de Lestang et al. 2012). Pueruli numbers are monitored at eight sites within the fishery every new moon (approximately monthly), with the lowest light and highest potential for collection (Phillips et al. 1978). These numbers form a puerulus index (PI) for fisheries management because, at this stage, they have ‘recruited’ into the western rock lobster fishery (de Lestang et al. 2012).

The role of physical oceanographic processes on the survival of P. cygnus phyllosoma during these vulnerable stages has been previously investigated through larval dispersal methods (Griffin et al. 2001; Feng et al. 2011; Caputi et al. 2018). Griffin et al. (2001) used Lagrangian methods in an individual-based model (IBM) in 2001. Modelling two sample years of low, high and average PI, they found a weak LC during early stages dispersed larvae west and north because of the wind. Conversely, as predicted pathways of puerulus crossed the continental shelf, a strong LC caused them to move south. Despite these differences and expected links to environmental factors, no changes were observed in the model PI (Caputi and Brown 1993; Caputi 2008). Therefore, the modelling by Griffin et al. (2001) suggested that biological responses from the hydrodynamic processes were the likely reason for the recorded natural variation in puerulus settlement. Subsequently, Feng et al. (2011) used six categories of larval stage P. cygnus to define these within an IBM. This model utilised Bluelink Reanalysis (BRAN) forecasting (Schiller et al. 2008) to assess varying oceanographic conditions as a form of larval dispersal for the larval-puerulus stage of P. cygnus. Caputi et al. (2018) built on this information by completing 24 model simulations to investigate the decrease in PI since the late 2000s. The model was forced with global and Australia-wide predicted atmospheric data with a resolution of 10-km grids (Caputi et al. 2018).

In 2008 and 2009, contrary to environmental conditions (Brown 2009), an unexpected decrease in the PI was observed, prompting a change in fishery management (de Lestang et al. 2015). After the PI recovered, different patterns in the puerulus settlement were observed. Differences include increased settlement in the north of the fishery (>31°S, Lancelin) and increased settlement in the second half of the season (November–March), as opposed to the previous peak being in October (Kolbusz et al. 2021a). A single cause for the unusually low puerulus settlement and its subsequent changes is unknown, and it is likely that multiple factors are involved. A period of sustained cooler water and low energetic conditions in the south-eastern Indian Ocean overlapped with the low settlement seasons (Kolbusz et al. 2021b). During early austral-summer months, there has been a consistent increase in the strength of the LC, coinciding with a decrease in northward wind-driven Cape Current (CC) c. 2008–09, which may explain the shifts in early settlement after this period (Pattiaratchi and Hetzel 2020; Pattiaratchi and Siji 2020; Kolbusz et al. 2021a). Whether these changes in settlement patterns have occurred primarily as a result of shifting oceanographic conditions or whether a biological mechanism, such as mature females spawning later, is yet to be fully understood.

Understanding why these changes have occurred for future climate-change scenarios and fishery management is essential. By using hindcast (in the past) oceanographic data, we can model the movement of particles and modify their movements to replicate those of a P. cygnus larvae at various developmental stages. As referred to previously, this has been completed using individual-based biophysical larval models (Feng et al. 2011; Caputi et al. 2018). Results showed several factors likely at play, one being onshore current systems mismatching with larvae leaving and returning to the coast (Caputi et al. 2018). Taking into account the finding of low energetic conditions in the south-eastern Indian Ocean (Kolbusz et al. 2021b), we have drawn our attention to the movement of the larvae with some biological considerations. We draw attention to the differences in larval movement before (2000–01), during (2007–09) and after (2010–13) the lower-than-expected puerulus settlement.

Our study aims to explain temporal variability (inter-annual and within a season) in settlement across latitudinal ranges for the early and late settlement periods of P. cygnus. We used three-dimensional (3-D) Lagrangian larval dispersal methods forced with 3-D advection fields from a numerical model. In particular, we simulated the dispersal of phyllosoma on hatching, transport in the open ocean, and pueruli movements towards recruitment areas. We also evaluated differences in transport temporally (early and late) and spatially (north and south) within each settlement season.

Off the western coast shelf of Australia, there is the southward-flowing Leeuwin Current (LC, present in the upper 0–300 m) and on the shelf is the northward-flowing wind-driven Capes Current (CC, surface) during the summer months (Fig. 2). The LC system includes the Leeuwin Undercurrent (LUC, >300 m) and the Ningaloo Current (NC) north of 24°S (Fig. 2; Woo and Pattiaratchi 2008; Pattiaratchi and Woo 2009). The LC is a poleward flowing eastern boundary current generated by an oceanic pressure gradient between north-western Australia and the Southern Ocean (Fig. 2c; Hamon 1965; Pearce 1991). This pressure gradient increases during a La Niña period, resulting in a stronger LC, warmer sea-surface temperature (SST) and increased rainfall over the Australian continent (Pattiaratchi and Buchan 1991; Feng et al. 2003; Telcik and Pattiaratchi 2014; Wijeratne et al. 2018; Pattiaratchi and Hetzel 2020; Pattiaratchi and Siji 2020). The opposite is true for El Niño. The mean southward volume transport of the LC peaks around 32.8°S, owing to inflow from the South Indian Counter Current (Wijeratne et al. 2018).

Fig. 2.  (a) A schematic diagram of the conditions along the western coast of Australia. Major currents and their approximate locations, which influence the transport of P. cygnus larvae, including the locations of key spawning areas in red. (b) Mean surface currents in austral summer, (c) mean surface currents in austral winter, (d) mean currents at 150-m depth in austral summer, (e) mean currents at 150-m depth in austral winter. Mean currents all calculated from ozROMS averaged velocities during 2000–18.

The LC, and associated currents, develop instabilities with the variations in topography over the extent of the current, forming eddies, meanders and offshoots (Andrews 1977; Pearce and Griffiths 1991; Cosoli et al. 2020). The region is associated with the most energetic eddy field along an eastern boundary (Pattiaratchi and Siji 2020). The Abrolhos Islands, the narrowing of the continental shelf slope south around 30°S, and the Perth Canyon off 32°S are major topographic features for the preferential generation of mesoscale eddies between 28 and 33°S (Fig. 2a; Feng et al. 2005; Rennie et al. 2007; Meuleners et al. 2008; Huang and Feng 2015; Cosoli et al. 2020). These eddies can persist for almost 8 months, with a mean radius of ∼100 km and generally retain their original formation (Fang and Morrow 2003; Cosoli et al. 2020). Contrary to other locations, the anticyclonic warm-core eddies of the south-eastern Indian Ocean have higher productivity than do cyclonic cold-core eddies because of entrainment of water from the continental shelf during their formation. As a result, eddies play a fundamental role in controlling species distribution and nutrient fluxes (Greenwood et al. 2007; Moore et al. 2007). The increasing strength of the LC is transferred to increases in meanders and eddies along the continental shelf (Feng et al. 2005).

Surface current variability on the continental shelf within the study region is predominantly wind-driven (Fig. 2a; Smith et al. 1991; Pattiaratchi et al. 1997; Rafiq et al. 2020). Southerly wind stresses overcome the alongshore pressure gradient during the summer months, moving upper layers offshore and favouring upwelling onto the continental shelf (Pearce and Pattiaratchi 1999). The CC is identified through these cooler waters, always inshore of the LC and usually inshore of the 50-m contour. It is initiated at approximately 34°S, flowing northward towards 27°S (Fig. 2a; Gersbach et al. 1999; Woo and Pattiaratchi 2008). Variations in the southerly wind stress and the Australasian monsoon season (January–March) over the austral summer cause offshore migration and weaker LC (Pearce and Pattiaratchi 1999; Woo and Pattiaratchi 2008; Feng et al. 2010; Pattiaratchi and Siji 2020). However, a weaker secondary peak in the LC can occur over December–January (Wijeratne et al. 2018). Thus, the LC migrates offshore over these predominantly sea-breeze-dominated summer months, whereas it floods the shelf during autumn–winter and dominates the distribution of water masses (Cresswell et al. 1989; Pattiaratchi et al. 1997; Woo and Pattiaratchi 2008).

Our Lagrangian larval dispersal model was forced with hourly advection fields from a 3-D hydrodynamic model on the basis of the Regional Ocean Modelling System (ROMS). ROMS (ozROMS, the ROMS configuration for Australia) was run in hindcast mode for Australia and the Indo-Pacific region between 2000 and 2018 (Wijeratne et al. 2018). OzROMS is a fully 3-D, high spatial-resolution (2–3 km) circulation model configured to include the entire Australian continental shelf, slope and the adjacent deep ocean by using ROMS. Further validation of the hindcast dataset is available in the extended validation section of the Supplementary material (Fig. S3–S11). Forcing for ozROMS includes tides, winds, air-sea exchanges and external forcing through the boundaries. These features are a distinct advantage for our application because prior larval dispersal studies of P. cygnus used a coarser resolution along the continental shelf. The details specific to the model set-up, hydrodynamic forcing and validation are found in Wijeratne et al. (2018).

Larval dispersal was completed using the open-source larval transport model Ichthyop (Lett et al. 2008). Ichthyop can track the trajectories of modelled larvae, or particles, in three dimensions through the ozROMS hydrodynamic model forcing at hourly intervals, including inputs for swimming particles, diurnal vertical migration and ontogenetic vertical migration (Lett et al. 2008; Robson et al. 2017). ‘Particles’, or in our case, simulated P. cygnus larvae, are released in the model domain, and the effects of the physical environment, particularly the currents, on the particles can be tracked.

Given the longevity of the P. cygnus larval phase and the uncertainty surrounding their ontogenetic vertical migration (OVM) patterns offshore, we broke down the larval dispersal model into the following three scenarios representative of the three life stages (Fig. 3): (1) hatching (Run H: early stage phyllosoma leaving shelf waters); (2) offshore (Run O: mid-stage phyllosoma); and, (3) settlement (Run S: early pueruli moving shoreward). We deliberately did not replicate the full larval phase because this has been addressed before (Griffin et al. 2001; Feng et al. 2011; Caputi et al. 2018) and is far more subjective to larval behavioural assumptions. Run H simulated newly hatched phyllosoma in their first 90 days, released on the sea bed at 3-km spacing between the 40 and 80-m depth contours, between 26 and 34°S, giving 58 520 particles released each season (Feng et al. 2011). Run O simulated phyllosoma already offshore for 60 days, released at 40-m depth at 120-km spacing in depths greater than 300 m over the whole model domain (22–36°S, out to 108°E), giving 20 070 particles released each season. Run S comprised the following two simulations: (1) puerulus once they have crossed the continental shelf (released on the 80-m contour) as passive particles; and, (2) puerulus at the continental shelf (released on the 200-m contour) as particles with directional swimming capabilities, giving 69 967 particles released each season. In both Run S simulations, particles were released at 40-m depth at a 3-km spacing between 26 and 34°S. The 3-km spacing was chosen on the basis of the size of the ozROMS horizontal resolution (Wijeratne et al. 2018).

Fig. 3.  Model scenario release and observation details for Run H (hatching scenario), Run O (offshore scenario) and Run S (settlement scenario).

Existing research has suggested that shifts in the oceanographic patterns over the past 20 years have contributed to the changes in puerulus settlement patterns and perhaps the unusually low settlement in 2008 and 2009 (Kolbusz et al. 2021a). Therefore, representative years for before (2000, 2001), during (2007–09) and after (2010–12) the settlement decrease, which intersected with ozROMS hourly data availability, were simulated for computation efficiency (Robson et al. 2017). The P. cygnus larval-cycle timing denotes that September to February of any given year is the hatching scenario data for the following settlement season (Fig. 3). Therefore, the settlement scenario is for the settlement season in question (Fig. 3). Season throughout refers to the puerulus settlement season (May–April).

All particles were forced by hourly ozROMS velocity fields in u and v directions and simulated using the Runge Kutta numerical advection scheme with a time step of 15 s. The small internal time step was chosen so that the Courant–Friedrichs–Lewy (CFL) criterion was not violated, causing instabilities when calculating w velocities. Hourly velocity fields were used, allowing for tides to be resolved, although the tidal range is small (∼0.5 m). The horizontal dispersion rate was set to 1 × 10⁻⁹ m² s⁻¹, as recommended within Ichthyop (Peliz et al. 2007; Lett et al. 2008). Particles were released every 10 days for the hatching (Run H) and offshore scenarios (Run O) and every 5 days for the settlement scenario (Run S; Fig. 3). The larval dispersal model was set to output every 12 and 24 h for Runs O and H respectively.

Within the settlement scenarios, particles were successfully ‘settled’ if they reached depths of less than 40 m, moved for at least 1 day, and stayed in depths of less than 40 m for at least 2 days. By contrast, if they had been moving for over 21 days and not crossed the 40-m contour, they were assumed to be perished because of their nektonic state and inability to survive (Feng et al. 2011).

Particles that moved outside the open ocean boundary (22–36°S, 108–117°E, Fig. 2) were considered out of bounds and removed. Particles were also removed if temperatures were outside survival bounds of 10–26°C. Particles that hit the coast were considered ‘beached’ and removed, unless it was within an appropriate settlement timeframe in Run S. If particles did not leave the continental shelf (>200 m) within 15 days of Run H, they were removed to account for high predation on the shelf (Caputi et al. 2018).

Ichthyop did not allow for both ontogenetic (OVM) and diurnal (DVM) vertical migration to co-occur. However, for the first 90 days of P. cygnus life, DVM occurs to a lesser extent, allowing only OVM methods to represent reality and testing against the passive Run H (Griffin et al. 2001; Feng et al. 2011). Ichthyop uses a method similar to that employed by the Connectivity Modelling System (CMS) to model OVM (Paris et al. 2017). The user inputs a probability matrix where the particles should be at time intervals after release. Following an approach similar to that of Feng et al. (2011), the particles stay between 0 and 20 m for the first 30 days. Between 30 and 90 days, the particles move between 0 and 80 m. Additional simulations were completed using OVM for Abrolhos and Big Bank (Fig. 2), key spawning areas (Feng et al. 2011). Particles were released in ozROMS at grid cells between 40- and 30-m depths at the Abrolhos Islands and between 100- and 80-m depth at Big Bank (Table 1).

Table 1.  Summary of performed model scenarios, purpose and associated sensitivity testing.

In Run S, swimming scenarios used the swimming action built within Ichthyop. This java code was modified to include a swimming direction instead of a random direction, which is the default. The direction remained at 64° from north towards the coast (Jeffs et al. 2005; Feng et al. 2011; Hinojosa et al. 2016). Three swimming speeds were tested, namely, 0.05, 0.1 and 0.15 m s⁻¹. Given that puerulus not likely to continuously swim at one speed (Phillips and Olsen 1975), 0.05 m s⁻¹ was considered the base Run S for subsequent follow-up results. As an additional sensitivity test, DVM was modelled for the Run S 0.05 m s⁻¹ scenario, with model puerulus near the surface at night and 20 m deep during the day (Feng et al. 2011).

The differences in the dispersal of particles after their first 3 months of life was the key response variable for Run H. Therefore, after 90 days, the locations of the particles for each iteration of Run H were collated into 1° bins. This was completed for early (September–December) and late (January–March) releases to correspond to the early and late spawning events, resulting in early and late settlements in the following year. The number of simulated phyllosoma within each degree bin was transformed into a percentage of the total present after 90 days to enable comparison among years.

For Run O, particle dispersal was mapped at 10, 30 and 60 days after release. These were transformed into a percentage of particles still present within the number of days to make all years comparable.

Each Run S output was processed by extracting only successfully recruited simulated puerulus that settled between 27.5 and 32.5°S to align with the locations of observed puerulus settlement (de Lestang et al. 2012). If particles settled after October, they were treated as settling in the latter portion of the season because of the changes in settlement patterns characterised by Kolbusz et al. (2021a). The latitudinal distribution of settled particles was explored through histograms for passive particles and those with swimming capabilities. These histograms were completed for the early and late portions of the season. Grouping the release latitude and settlement latitudes of model puerulus into connectivity matrices gave a spatially comparable indication of how the particles moved before successful settlement.

Given the results of the connectivity matrices, the shift between the release and settlement latitudes was explored further. Further exploration was completed only for Run S 0.05 m s⁻¹ scenario because the number of successful particles settling was higher than passive, but less different when the speed was increased. The following three fixed factors were used to test changes in the mean latitude shift from release to settlement location before, during and after the low puerulus settlement: in the early or late season (two levels; early and late), the location of settlement (two levels; north, >31°S; and south, <31°S) and before, after or during the period of low settlement (three levels; before, during and after). For example, if particles settled north of 31°S, they were said to have settled north of the fishery (Kolbusz et al. 2021a). The mean latitude shift for each successful model puerulus recruit was averaged for both the early and late settlement season and separated into the northern and southern settlement. The years were treated as random and used as replicates nested within period (eight levels; 2000–01, 2007–12). An ANOVA was conducted to test for changes in the mean latitude shift of particles using the PERMANOVA+ add-on package for PRIMER (ver. 6, see https://www.primer-e.com/, using 9999 permutations; Anderson et al. 2008) and a Euclidean distance measure. Significant terms and interactions were investigated using a posteriori pairwise comparison (see Anderson 2001).

A comparison of DVM particles was completed with the Run S 0.05 m s⁻¹ as the baseline. In addition, the contrast between the DVM and non-DVM particles was explored by dividing them into early or late successful settlements, as was completed for all Run S simulations.

Only ‘alive’ particles, defined as those currently active, were considered in the analysis for all scenarios. There was an average of 88% mortality within Run H in passive-transport simulations and 85% in the OVM transport simulations (Supplementary Table S1). There was an average of 44% mortality within Run O (Supplementary Table S2) and 70% mortality or not-recruited within Run S (Supplementary Table S3).

Run H found that during and after the period of unusually low settlement, particles were dispersing to a greater extent in both early and late releases (Fig. 4). However, in 2001 and 2007, the higher-density area was further north (Fig. 4). After 2007, the late release particles shifted southward. Before the low puerulus settlement, there were specific areas of higher-density particle retention, which were likely to be within eddy features (Fig. 4). Both passive and ontogenetic verticle migration (OVM) particles had high mortality for beaching in 2007 and 2008 (38%; Table S1).

Fig. 4.  Distribution of particles (assumed phyllosoma) after 90 days of passive transport from hatching areas (Run H). White values indicate the number of active particles used for the associated timeframe. Before the decline simulations: (a, b) 2000 releases, 2000 releases (c, d) 2000 releases. During the decline simulations: (e, f) 2007 releases; (g, h) 2008 releases. After the decline simulations: (i, j) 2011 releases; (k, l) 2012 releases. Early releases are from September to December of a given year, late releases are January–February of the year after.

After 90 days, the locations of OVM particles were reasonably similar to those transported passively, with both displaying less concentration of particles after 2007 (Fig. 5g–l). After this time, dispersal to the north and south is more extensive for early and late releases throughout the unusually low settlement. There was a considerable dispersal of particles to the north in 2007, dissimilar to the corresponding passively transported simulation. Interestingly, in 2012 after the period of unusually low settlement, the higher-density locations were close to the shelf. There is a higher-density offshore at 33−34°S, but it is still more dissipated than earlier.

Fig. 5.  Distribution of particles (assumed phyllosoma) after 90 days of OVM transport from hatching areas (Run H). White values indicate the number of active particles used for the associated timeframe. Before the decline simulations: (a, b) 2000 releases; (c, d) 2001 releases. During the decline simulations: (e, f) 2007 releases; (g, h) 2008 releases. After the decline simulations: (i, j) 2011 releases; (k, l) 2012 releases. Early releases are from September to December of a given year, late releases are from January to February of the year after.

Run H particles released from the Abrolhos and Big Bank with OVM show similar dispersal patterns in the early releases (Fig. 6), aside from 2012, where wider dispersal occurred. In the early releases of 2007, a large portion (∼10%) of particles ended up at 32°S offshore, which is inconsistent with the other years (Fig. 6e). Late releases (Fig. 6) showed higher spatial variation than did early releases, with the highest-density area moving further north and inshore with each simulation. Still, by 2012, the particles were widely dissipated.

Fig. 6.  Distribution of particles (assumed phyllosoma) after 90 days of larval dispersal with OVM transport from Abrolhos (ABR) and Big Bank (BB), which are key hatching areas (Run H). White values indicate the number of active particles used for the associated timeframe and figure. Before the decline simulations: (a, b) 2000 releases; (c, d) 2001 releases. During the decline simulations: (e, f) 2007 releases; (g, h) 2008 releases. After the decline simulations: (i, j) 2011 releases; (k, l) 2012 releases. Early releases are from September to December of a given year, late releases are from January to February of the year after.

Offshore releases of particles, representing the mid-stage phyllosoma life-stage over the winter months, showed only minor differences among years (Fig. S1). Despite the same number of particles, a large portion (∼20%) left the boundary in 2000 after only 10 days. Approximately 5% of particles died as a result of being outside of the temperature bounds (Table S2). In 2011, a higher particle density was along the shelf region at 34°S after 30 days. Although by 60 days, the density of particles increased in areas off the shelf up to 112°E (Supplementary Fig. S1). The same patterns were seen in 2000. There was a reduction in relative density in the northern-most sections during 2011 and 2012 compared with the earlier years (Fig. S1).

The maximum number of particles to settle in Run S (passive) was in 2009 (4362, Fig. 7a), with the majority of particles settling in the earlier portion of the season. However, there was an increase in settlement in the latter portion of the season after 2001 (Fig. 7a). An average of 74% of particles did not reach the nearshore, with 20% outside of temperature bounds (Table S3).

Fig. 7.  Proportion of particles released that settled within the early and late portions of the seasons (Run S). (a) Passive particles. (b) Swimming particles at 0.05 m s⁻¹. (c) Swimming particles at 0.01 m s⁻¹. (d) Swimming particles at 0.15 m s⁻¹. The number above the bars are the total number of particles settled for the simulation (early and late). Early settlement is a successful settlement between September and October. Late settlement is a successful settlement between November and March.

Before and after the low puerulus settlement, the peak latitude of settlement distribution peaked at 28.5°S. During the low settlement, this peak was further south at 30°S (Fig. 8a). The latitudinal distribution of the late settlement was similar to that of the early settlement (Supplementary Fig. S2). In both early and late portions of the seasons, there was also a drop in the settlement between 29.5 and 29°S each year (Fig. 9a, Fig. S2).

Fig. 8.  Settlement latitudes for successfully settled particles in the early portion of the season (Run S). Early settlement is a successful settlement between September and October. (a) Passively transported particles. (b) Swimming particles at 0.05 m s⁻¹. (c) Swimming particles at 0.1 m s⁻¹. (d) Swimming particles at 0.15 m s⁻¹.

Fig. 9.  Connectivity matrices (as a percentage) of successfully settled particles in the early and late portions of the season (Run S). Early settlement is a successful settlement between September and October. Late settlement is a successful settlement between November and March. The black line corresponds to where the source (release) and sink (settlement) latitudes are the same. (a) Passively transported particles. (b) Swimming particles at 0.05 m s⁻¹.

The number of simulated puerulus settling increased when swimming capabilities were added (Fig. 7b–d), with the conditions in the early portion of the season (September–November) being more favourable than those later (December–March) in the season (Fig. 8c–d). However, an increase in swimming speed did not considerably increase the survival of particles (Fig. 7, Table S3). Furthermore, settlement decreased over the early and late portions of the season in 2007 and 2008 and was most prominent for the 0.05 m s⁻¹ simulations (Fig. 7b).

The latitudinal distribution of modelled settlement occurred in two prominent peaks, similar to the non-swimming scenario (Fig. 8a). However, the inclusion of swimming caused the distributions to widen into a more uniform spread across the latitudes. This change was apparent as the swimming speed increased. In most scenarios and seasons (including passive), very little settlement occurred north of 28°S, and there was also less settlement just south of 29°S. The latter portion of the season had similar distribution patterns for all swimming particles (Fig. S2).

For the passive model puerulus, during the early portion of the season, the most important release region was ∼27 to 28°S (Fig. 9a). However, this was further south for the latter portion of the season, ∼30 to 31°S for 2000–01 and 2007–09. It shifted to be more similar to the early portion with a northern release latitude (Fig. 9a). The spread of the release-settlement groupings moves further below the black line (where release and settlement latitudes are the same) over time in both the early and later portions of the season (Fig. 9a). As a result, this shows that the particles are being transported further South before eventually settling.

When swimming speed was added (0.05 m s⁻¹), as particles moved directly toward the coast, their release and settlement latitudes were more similar (black line, Fig. 9b). The most important release latitude was 30–32°S for early settlement. This was less distinguishable for the latter settlement, where the latter settlement grouping was more dispersed than the earlier (Fig. 9b). However, both early and late portions of the season showed that the shift in latitude between release and settlement moved southward during and after the period of observed reduced settlement (Fig. 9b).

A comparison of the latitudinal shift between the release and settlement latitude was calculated for all 0.05 m s⁻¹ swimming particles (Supplementary Table S4). Before the low settlement, particles from the north and south were typically transported northward (Fig. 10a). However, during and after the unusually low puerulus settlement, particles were typically transported southward if they originated in the north (during t = 5.87, P_mc < 0.01, after t = 3.999, P_mc < 0.05), and more weakly northward if they originated in the south (during t = 8.153, P_mc < 0.005, after t = 7.024, P_mc < 0.01; Fig. 10a). Over the whole simulation period, transport in the north was not significantly different over the early and late portions of the season. However, in the south, the source–sink pattern of the late portion of the season changed to a significantly northward transport (t = 4.512, P < 0.005; Fig. 10b).

Fig. 10.  Mean latitude shift, in degrees, between the release and settlement latitudes for particles swimming at 0.05 m s⁻¹ (Run S) (a) for before (2000–01), during (2007–09) and after (2010–12) the decline in the north (>31°S) and south (<31°S). Negative is, therefore, a movement south. The interaction between location (north or south) and period (before, during or after) was significant in the ANOVA model (b) for the early (September–October) and late (November–March) portions of the season in the north (>31°S) and south (<31°S). The interaction between location (north or south) and timing (early or late) was significant in the ANOVA model. Results of pair-wise comparisons are indicated by lower-case letters in a and b, with different letters indicating a significant difference. Positive is a northward shift and negative is a southward shift.

Overall, particles were transported more southwards in the north during and after the period of low settlement (Fig. 10). Trajectory mapping further supports this pattern, where settled particles in a predominant location in the north are more likely to originate from further northwards during and after the period of low settlement (Fig. 11). In comparison to the south, where before the period of low settlement, particles typically originate from a range of latitudes along the continental shelf in both early and late releases (Fig. 12); however, over time there has been a tendency for particles to originate from further north, particularly in the late releases (Fig. 12).

Fig. 11.  Drift trajectories for a subset of successfully settled particles that settled in the northern peak latitude area between 28.6°S and 29.1°S (Run S). Blue is their start location and red is their end location. Red values indicate the number of successfully settled particles used for the associated timeframe and figure. Before the decline simulations: (a, b) 2000 releases; (c, d) 2001 releases. During the decline simulations: (e, f) 2007 releases; (g, h) 2008 releases; (i, j) 2009 releases. After the decline simulations: (k, l) 2010 releases; (m, n) 2011 releases; (o, p) 2012 releases. Early settlement is from September–October of a given year, late settlement is November–March of the year after.

Fig. 12.  Drift trajectories for a subset of successfully settled particles that settled in the southern peak latitude area between 31.5 and 31.9°S (Run S). Blue is their start location and red is their end location. Red values indicate the number of successfully settled particles used for the associated timeframe and figure. Before the decline simulations: (a, b) 2000 releases; (c, d) 2001 releases. During the decline simulations: (e, f) 2007 releases; (g, h) 2008 releases; and (i, j) 2009 releases. After the decline simulations: (k, l) 2010 releases; (m, n) 2011 releases; (o, p) 2012 releases. Early settlement is from September–October of given year, late settlement is November–March of the year after.

The pattern of settlement, allowing for diurnal vertical migration (DVM) behaviour of particles (Fig. 13), was relatively comparable to simulations without DVM (Fig. 7b). However, on average, there was 15% less settlement in later periods than earlier in the settlement season (Fig. 13). In addition, the late settlement showed a clear declining pattern until 2007 and then incrementally recovered to past the previous maximum by 2012 (Fig. 13). Again, this varies from the simulations without DVM (Fig. 7b), where the late settlement exhibited as much temporal variation. An additional difference was an abrupt decline in the early settlement in the 2010 season (Fig. 13), which was not as pronounced without DVM (Fig. 7b).

Fig. 13.  Proportion of particles released that settled within the early and late portions of the seasons (Run S, as in Fig. 7) for swimming particles at 0.05 m s⁻¹ with diurnal vertical migration. The numbers above the bars are the total number of particles settled for the simulation (early and late). Early settlement is a successful settlement between September and October. Late settlement is a successful settlement between November and March.

Our results suggest that the variations in observed puerulus settlement between 2000 and 2012 are a product of variation in the localised driving forces of the Leeuwin (LC) and Capes (CC) currents over the summer months. The years of low settlement coincided with a low energy period in the south-eastern Indian Ocean, where there were neutral ENSO conditions over an extended period (Pattiaratchi and Siji 2020; Kolbusz et al. 2021b). In 2011, when there was a strong La Niña, recovery in puerulus numbers began, but not to the expected extent based on the strength of the LC (de Lestang et al. 2015). This is surmised in Fig. 14. We have demonstrated that the cause of settlement failure could be identified by isolating the different early life-cycle phases of P. cygnus and the forcing mechanisms at each dispersal phase.

Fig. 14.  Schematic summarising the changes in the movement of puerulus onshore before, during and after the period of low settlement. Including the change in the wind-driven Capes Current strength influenced by the increase in the Leeuwin Current. Arrows are for visual representation are purposefully approximate.

Our settlement results (Run S) have shown that particles are being transported from locations substantially and significantly further north during the period of unusually low settlement compared with prior. This trend also continued after the period of unusually low settlement (2010–12). The latitudinal shift for north-settling particles was similar in the early and late portions of the season. However, particles moved significantly further north to settle during the latter part of the season in the south. This matches observed changes in puerulus settlement, with a noticeable drop in the PI in the southern sites during and after the recruitment failure (Kolbusz et al. 2021a). In 2010 (the late portion being in the early 2011 summer) and 2011, the LC was strong over the later portions of those seasons, thereby weakening the CC by flooding the shelf (Benthuysen et al. 2014); this is likely to have decreased the number of observed puerulus compared with what was typically expected during a strong La Niña (Caputi 2008). The difference in the early and late latitudinal shifts also matched the differences in recovery of the puerulus index after the period of unusually low settlement experienced between the north and south of the fishery (Kolbusz et al. 2021a).

Whether north or south, an increase in the degree difference between the release and settlement latitude may imply that the distance travelled by the model pueruli has increased. The distance travelled for a puerulus is highly important for its survival because of its nektonic state (Lemmens 1994; Fitzgibbon et al. 2014). The further they travel, the higher the utilisation of energy reserves. With the increased dispersal of particles shown in our hatching runs, phyllosomas reaching the continental shelf to undergo metamorphosis may already have depleted energy reserves if they were not within eddy systems (Saunders et al. 2012; Waite et al. 2019). Model pueruli that successfully settled in the south of the fishery, particularly during the late portion of the season, travelled further to reach settlement locations. In the north, despite there being a southward shift in the direction they travelled to get to the settlement location (because of the LC summer increase), the distance travelled was more uniform over the periods. The differences in distance travelled may explain the more significant recovery in PI observed in the northern parts of the fishery (Kolbusz et al. 2021a).

The division of the larval cycle into three simulations (Run H, Run O and Run S), drew attention to the dispersal dynamics of P. cygnus larvae, bypassing biological considerations such as temperature-dependent growth (Liddy et al. 2004; Steele and Neuheimer 2022), primary productivity (Caputi et al. 2001) and food availability (Waite et al. 2019). On the eastern coast of Australia, decreases in the pelagic larval duration for Melicertus plebejus (prawn; Everett et al. 2017) and Pomatomus saltatrix (bluefish; Schilling et al. 2020) have been due to degree-day growth changes, in combination with mesoscale processes. Complete biophysical models considering degree-day growth have also been completed for P. cgynus (Feng et al. 2011; Caputi et al. 2018); however, the culmination of this research and the above considerations highlight the complexities surrounding the years of low settlement of P. cgynus. The difference in the movement of particles during and after the period of low settlement, combined with an increase in temperature from the LC, would have reduced the intermoult period and possibly their size (Liddy et al. 2004). Therefore, a decrease in their pelagic larval duration during early settlement could have led to a decreased CC (pre-November), possibly discouraging the larvae off the shelf.

After 3 months, the wider dispersal of particles by 2012 highlights that there has been less entrainment of particles within eddy structures, likely reducing their return pathway to the coast. The entrainment of phyllosoma and other long-lived larvae within eddies is regarded as positive for prey availability and transport (Mullaney and Suthers 2013; Cetina-Heredia et al. 2019; Woodings et al. 2019; Clavel-Henry et al. 2020). It has previously been suggested that the period of unusually low settlement may have occurred because of a ‘rogue’ eddy system entraining newly hatched larvae from the key spawning areas and transporting them offshore (Brown 2009; Säwström et al. 2014; Waite et al. 2019). However, Run H has shown that this is unlikely to have occurred. Instead, there is a tendency for particles to become increasingly dispersed leading up to, during and after the period of low settlement, with limited evidence of entrainment in a single eddy structure.

Particles released offshore during austral winter (Run O) were also influenced by the strength of the LC and an associated increase in eddy kinetic energy (EKE). The number of particles lost from the simulations was highest in 2000 when EKE in the system peaked (Pattiaratchi and Siji 2020; Kolbusz et al. 2021b). There was a similarly increasingly southern distribution of particles in 2011 to 2000, which coincided with another peak in EKE (Pattiaratchi and Siji 2020; Kolbusz et al. 2021b). The only simulation where we observed an increased density of particles overlapping with the flow of the LC was after 30 days in 2011 when the LC was undoubtedly the strongest (Boening et al. 2012; Benthuysen et al. 2014; Wijeratne et al. 2018). Alongside the increased dispersal in Run H, the strong LC in 2011 may have been too strong for phyllosoma to survive. Fisheries managers expected a significant increase in settlement given the strong La Niña conditions in 2011 (de Lestang et al. 2015). However, only a small increase in PI was observed (Kolbusz et al. 2021a). The LC also had a more considerable spatial extent (Huang and Feng 2015) throughout the unusually low settlement, which may have resulted in inadequate energy in the system to increase the PI substantially.

Despite minor inter-annual differences, the number of particles that successfully settled during Run S remained relatively consistent. This was expected because we did not model transport over the complete life-cycle of the larvae (but see Feng et al. 2011; Caputi et al. 2018). A slight decrease was recorded over 2007 and 2008, in line with the period of low settlement when swimming was added (Fig. 7b); however, not in 2009, where a low PI was also recorded (de Lestang et al. 2015; Kolbusz et al. 2021a). The lower number of successfully settled particles in the late portion of the season in Run S, including the diurnal vertical migration (DVM) sensitivity runs, suggests that the weak secondary LC peak over January and February inhibits puerulus from reaching the nearshore. The same is not likely to occur over the early portion of the season where the CC is flowing northward, and LC would usually be weaker (Gersbach et al. 1999; Wijeratne et al. 2018).

Separating the successfully settled particles from Run S into their locations along the coast showed areas where shelf and nearshore processes assist puerulus transport to likely settlement locations. More temporal variability in settlement latitudes was seen within the early settlement than the late, although, overall, they showed similar regions of increased settlement. Both passive- and swimming-particle simulations indicated a region in the north (28.2–28.8°S), where settlement decreases throughout the low settlement period, causing the distribution to skew southward. This north region includes Dongara, where the PI has historically been the highest and decreased significantly during the period of unusually low settlement (de Lestang et al. 2015). This suggests that the transport conditions during the period of the unusually low settlement were not conducive to transporting puerulus towards these key latitudes within the fishery.

Consistent with the literature on other spiny lobster species, our research showed the influences of nearshore oceanographic processes on the settlement of pueruli. For Jasus edwardsii, south-westerly wind-driven surface currents also assist transport to the South Australian coastline (Hinojosa et al. 2017). In addition, hurricanes and their associated sea-level changes have been correlated with Caribbean spiny lobster (Panurlius argus) settlement patterns, suggesting that the conditions assist the movement of recruits onshore (Briones-Fourzán et al. 2008). Throughout the unusually low settlement for P. cygnus, the boundary current, the LC, was stronger over the summer, flooding the shelf, and reducing the CC strength. At the same time, it functions to assist puerulus settlement. The same occurs within the East Australia Current, where the jet-like stream can act as a barrier to the settlement of Sagmariasus verreauxi (Roughan et al. 2011; Cetina-Heredia et al. 2015). In essence, it is a combination of offshore abundance and condition, followed by active swimming and oceanographic shelf processes that drive the successful settlement of P. cygnus and other spiny lobster species globally (Yeung et al. 2001; Goldstein and Butler 2009; Linnane et al. 2010; Singh et al. 2018).

An increase in the dispersal properties of the LC during hatching, alongside its increased strength in summer, while the CC decreased, is likely to have given rise to alternative pathways for puerulus reaching the nearshore after the low settlement period. It is important to note that with only minimal field observations, little is definitively known about the movements of P. cygnus during its larval duration, in particular, cues for metamorphosis into puerulus and likely swimming speeds (Jeffs et al. 2005; Phillips and McWilliam 2009). This is a limitation of our study, alongside the assumptions of larval dispersal modelling. In particular, biological processes such as variation in productivity and predation were not factors in this analysis aside from a timing constraint (40 days) to leave the shelf in Run H. However, we aimed to highlight the possible influence that oceanographic conditions could have on planktonic and weakly swimming particles before, during and after this period of unusually low settlement.

Our results have added to the growing body of knowledge surrounding the influence of the LC and coastal currents on P. cygnus larvae and settlement patterns represented by PI (Caputi and Brown 1993; Pearce and Phillips 1988; Kolbusz et al. 2021b). Our research indicated that localised influences of the LC interacting with CC can positively and negatively affect the successful settlement of P. cygnus puerulus. However, this depends on where and when it has reached the shelf and whether it is in the early or late portion of the season. This research, therefore, calls for greater monitoring of the CC strength when considering larval stages of species over the summer months along the western coast of Australia, for example, characterisation of the spatial scale of the LC and CC from satellite images in austral winter and summer. Importantly, larval dispersal modelling should consider it necessary for the CC to be resolved.

Our findings have significant implications on fishery management and understanding of the transport mechanism of the early life of the western rock lobster. The southward-flowing Leeuwin Current (LC) and the northward-flowing Capes Current (CC) on the shelf interact over summer to generate different pathways for puerulus to reach nearshore areas. Before the unusually low puerulus settlement period, the LC was weaker, therefore not impeding the CC northward shelf flow. This allowed pueruli to be transported more directly to suitable habitats. During and after the period of low settlement, puerulus reached nearshore areas more successfully from further north, particularly in the southern latitudes of the fishery. The same occurred afterwards, although to a lesser extent. These differences were highlighted statistically and through trajectory mapping. We have shown that it is not so much the strength of the LC but somewhat nuanced pathways influencing the CC in summer, that characterises puerulus settlement into the fishery.

Supplementary material is available online.

Monthly, daily and hourly mean 3-D currents, temperature and salinity together with mean water levels available for the ozROMS simulations described in this paper are available through contacting the first author.

The authors declare that they have no conflicts of interest.

J. Kolbusz was supported by the UWA University Postgraduate Award and Australia Research Training Program scholarships.