DOI: 10.5061/dryad.2547d7x23
title: Using predictive models to identify kelp refuges in marine protected areas for management prioritisation

description: Kelp forests serve as the foundation for shallow marine ecosystems in many
 temperate areas of the world but are under threat from various stressors,
 including climate change. To better manage these ecosystems now and into
 the future, understanding the impacts of climate change and identifying
 potential refuges will help to prioritise management actions. In this
 study, we use a long-term dataset of observations of kelp percentage cover
 for two dominant canopy-forming species off the coast of Victoria,
 Australia: Ecklonia radiata and Phyllospora comosa. These observations
 were collected across three scuba sampling programs that extend from 1998
 to 2019. We then associated those observations with habitat and
 environmental variables including depth, seafloor structure, wave climate,
 currents, temperature, and population connectivity in generalised additive
 mixed effects models and used these models to develop predictive maps of
 kelp cover across the Victorian MPAs. These models were also used to
 project kelp coverage into the future by replacing wave climate and
 temperature with future projections (2090, RCP4.5 and RCP8.5). Once all
 the past, current, and future predictions were compiled, we calculated
 percent cover change from 1998-2019, stability over the same period, and
 future predicted change in percent cover (2019-2090) to understand the
 dynamics for each species across the MPAs. We also used the current
 percentage cover, stability, and future percentage cover to develop a
 ranking system for classifying the maps into high climate vulnerability,
 climate vulnerability, neutral, potential refugia, and likely refugia. A
 management framework was then developed to use those refugia ranking
 values to inform management actions and we applied this framework across
 three case studies: one at the scale of the MPA network and two at the
 scale of individual MPAs, one where management decisions were the same for
 both species and one where the actions were species-specific. This study
 shows how species distribution models, both contemporary and with future
 projections, can help to identify potential refugia areas that can be used
 prioritise management decisions and future-proof restoration actions.
description: <strong>Parks Victoria Subtidal Reef Monitoring Program
 (SRMP)</strong> Parks Victoria SRMP targets areas
 within MPAs and sites directly adjacent to them and uses visual census
 methods developed by Edgar and Barrett (1997), Edgar and Barrett (1999),
 and Edgar et al. (1997). There are 80 SRMP sites used in this study that
 extend from 1998 to 2013. Each site within the sampling program was marked
 by GPS and visited in subsequent survey years where a 200 m transect was
 laid along a shallow (&lt; 10 m) contour. To quantify the percentage
 coverage of kelp, 0.25 m<sup>2</sup> quadrats, which are
 divided into a grid of 7 x 7 perpendicular wires to provide 49 points plus
 one corner making 50, are placed at 10 m intervals along the transect. The
 cover of kelp is then estimated by counting the number of times the
 species intersects with one of the points. These point counts are then
 divided by 50 and averaged across all quadrats in the transect to come up
 with a percentage coverage of each kelp for each site. For this project,
 only the percentage cover of <em>E. radiata</em> and
 <em>P. comosa</em> were extracted from the dataset.
 <strong>Reef Life Survey (RLS)</strong>
 Reef Life Survey (RLS) is a citizen science program where
 volunteer SCUBA divers conduct visual surveys of underwater reefs. This
 program began surveying 95 sites along the Victorian coast in 2007 and
 continues to extend the subtidal biological observation time series,
 including to 2019, the last year used in this study. See Edgar and
 Stuart-Smith (2009) for detailed methods. Briefly, a GPS point is taken at
 each site and it is surveyed using a 50 m line transect along a shallow
 depth contour with two or more contours targeted at each site. To collect
 data on kelp coverages, digital photo quadrats taken downward from a 50 cm
 height above the seabed are acquired at 2.5 m intervals along the transect
 line. Post-processing of these photo quadrats requires classifying the
 images into 16 functional or morphological categories of algae using
 CATAMI (Althaus et al. 2015), which is the standard image classification
 system for Australia. Once classified, the point counts of <em>E.
 radiata</em> and <em>P. comosa</em> were used to
 calculate percentage cover of each species.
 <strong>Victorian Fisheries Authority
 (VFA)</strong> To increase the coverage and time
 series of kelp observations, the VFA data were used. For this dataset,
 diver surveys are conducted annually at 202 monitoring sites across the
 coast of Victoria. These surveys began in 2002 and continued until 2016.
 At each site, the divers estimate the percent cover of <em>E.
 radiata</em> and <em>P. comosa</em> within six 30 m by 1
 m transects at random cardinal directions from central site coordinates
 taken using a GPS point. The estimates for percent cover are based on
 visual approximation within 30 sections along each transect and then
 averaged to the nearest 10% across all transects for each site.
 Data from these three sampling programs were collated into a
 single dataset for each species: <em>Ecklonia radiata</em> and
 <em>Phyllospora comosa </em>for statistical analyses. The
 resulting compiled dataset included over 300 locations sampled over 22
 years across Victoria, both inside and outside marine protected areas.
 Each sampling program captured data using SCUBA following slightly
 different sampling protocols (see Supplementary Information 1.1 for
 details on differences in sampling). Then, for the sampling methods
 (either using quadrats, photo quadrats, or multiple transects), an average
 percentage coverage for each species was calculated per site per year
 (i.e., all quadrats or transects at each site were averaged to produce one
 percentage cover value per site), resulting in a time series of percentage
 cover at each site for both species.
description: # Using predictive models to identify kelp refuges in marine protected
 areas for management prioritisation
 [https://doi.org/10.5061/dryad.2547d7x23](https://doi.org/10.5061/dryad.2547d7x23) ## Description of the data and file structure ### Percentage cover of Ecklonia radiata and Phyllospora comosa across Victoria, Australia and associated environmental covariates Access this dataset on Dryad ([https://doi.org/10.5061/dryad.2547d7x23](https://doi.org/10.5061/dryad.2547d7x23)) These datasets contain separate percentage cover information for *Ecklonia radiata* and *Phyllospora comosa* along the coastal waters of Victoria, Australia along with the habitat and environmental conditions associated with the spatial location of each sample. This dataset is a combination of three different sampling programs, outlined in the **Methods** section. ## Description of the data and file structure These datasets are .csv files that contain information on percentage cover of *Ecklonia radiata* (Ecklonia_radiata.csv) and *Phyllospora comosa* (Phyllospora_comosa.csv) at locations across Victoria, Australia. The location information is stored both as Latitude and Longitude (WGS84) and Eastings and Northings (GDA1994 VICGRID). These datasets were used to run the gams models with the percentage cover as the response variable and several of the covariates as explanatory variables. The covariates were extracted from raster datasets representing each of them. A list of the column names and their meanings are on the next page. ## Sharing/Access information You can access the data here: [https://portal.aodn.org.au/](https://portal.aodn.org.au/) You can the bathymetry data for deriving the physical variables: [https://portal.ga.gov.au/persona/marine](https://portal.ga.gov.au/persona/marine) ## Code/Software Example gam model run for *Ecklonia radiata* using the one of the datasets and the mgcv package in R: Erad.n250 &lt;- gam(mean_PercentCover ~ s(s_sst_c) + s(depth, k=2) + s(w_mov, k = 2) + s(InDegree, k=2) + fBioUnit_I + s(slopeposition_n250, k=5) + s(landform_n250, k=5) + s(a_acs, k=5) + s(roughness_n250, k=3) + s(LocalRetention) + s(landform_n250, k=3) + s(surfacereliefratio_n250) + s(fYear, bs = 're'), family=betar(link = "logit"), data = Erad_train, method = 'REML') ### Files and variables #### File: Ecklonia\_radiata.csv **Description:**  ##### Variables * OID_: object ID - Location: Bioregional unit the sampling occurred in - Species: The species name of the macroalgae - SiteName: The survey site name - Year: the year the data were collected - SiteCode: Unique site code for each survey site - Latitude: Latitudinal location of the site - Longitude: Longitudinal location of the site - n: number of samples - mean_PercentCover: Average percentage cover across samples - meanPercentCoverI: Percent cover converted to an integer value - BioUnit: The bioregional unit associated with the site location - MPAorRef: Either the name of the MPA the site was collected in or a "Ref" area associated with an MPA - depth: Depth of the site (m) - surfacereliefratio_n3: measure of surface relief using GA-SaMMT v1.2 with a neighbourhood of 3 cells (no units) - slopeposition_n10: Relative elevation of a cell compared to the overall seascape calculated using the benthic terrain modeller with a neighbourhood of 10 cells (no units) - slopeposition_n20: Relative elevation of a cell compared to the overall seascape calculated using the benthic terrain modeller with a neighbourhood of 20 cells (no units) - slopeposition_n100: Relative elevation of a cell compared to the overall seascape calculated using the benthic terrain modeller with a neighbourhood of 100 cells (no units) - surfacereliefratio_n10: measure of surface relief using GA-SaMMT v1.2 with a neighbourhood of 10 cells (no units) - surfacereliefratio_n20: measure of surface relief using GA-SaMMT v1.2 with a neighbourhood of 20 cells (no units) - surfacereliefratio_n250: measure of surface relief using GA-SaMMT v1.2 with a neighbourhood of 250 cells (no units) - slopeposition_n500: Relative elevation of a cell compared to the overall seascape calculated using the benthic terrain modeller with a neighbourhood of 500 cells (no units) - surfacereliefratio_n500: measure of surface relief using GA-SaMMT v1.2 with a neighbourhood of 500 cells (no units) - landform_n3: Morphology of the seafloor based on the size, slope, length-width ratio using GA-SaMMT v1.2 with a neighbourhood of 3 cells (no units) - landform_n100: Morphology of the seafloor based on the size, slope, length-width ratio using GA-SaMMT v1.2 with a neighbourhood of 100 cells (no units) - landform_n250: Morphology of the seafloor based on the size, slope, length-width ratio using GA-SaMMT v1.2 with a neighbourhood of 250 cells (no units) - landform_n500: Morphology of the seafloor based on the size, slope, length-width ratio using GA-SaMMT v1.2 with a neighbourhood of 500 cells (no units) - landform_n10: Morphology of the seafloor based on the size, slope, length-width ratio using GA-SaMMT v1.2 with a neighbourhood of 10 cells (no units) - landform_n50: Morphology of the seafloor based on the size, slope, length-width ratio using GA-SaMMT v1.2 with a neighbourhood of 50 cells (no units) - roughness_n10: A measure of seafloor complexity calculated using the variation in slope and aspect vectors calculated using Benthic Terrain Modeller with a neighbourhood of 10 cells (no units) - roughness_n250: A measure of seafloor complexity calculated using the variation in slope and aspect vectors calculated using Benthic Terrain Modeller with a neighbourhood of 250 cells (no units) - roughness_n50: A measure of seafloor complexity calculated using the variation in slope and aspect vectors calculated using Benthic Terrain Modeller with a neighbourhood of 50 cells (no units) - roughness_n100: A measure of seafloor complexity calculated using the variation in slope and aspect vectors calculated using Benthic Terrain Modeller with a neighbourhood of 100 cells (no units) - roughness_n20: A measure of seafloor complexity calculated using the variation in slope and aspect vectors calculated using Benthic Terrain Modeller with a neighbourhood of 20 cells (no units) - s_sst_c: Average summer sea surface temperature. SST values are associated with the year of survey (degrees Celsius) - a_acs: Average annual current speed. Current speed values are associated with the year of survey (m/s) - w_mov: Average winter wave orbital velocities. Orbital velocity values are associated with the year of survey (m/s) - LocalRetention: The proportion of released spores that settle back to the source site from the biophysical modelling. Values are associated with the year of survey (no units) - InDegree: The total number of significant connections coming into a site from the biophysical modelling. Values are associated with the year of survey (no units) - Easting_VG94: Easting value from the VICGRID 1994 coordinate system (m) - Northing_VG94: Northing value from the VICGRID 1994 coordinate system (m) #### File: Phyllospora\_comosa.csv **Description:**  ##### Variables * OID_: object ID * Location: Bioregional unit the sampling occurred in * Species: The species name of the macroalgae * SiteName: The survey site name * Year: the year the data were collected * SiteCode: Unique site code for each survey site * Latitude: Latitudinal location of the site * Longitude: Longitudinal location of the site * n: number of samples * mean_PercentCover: Average percentage cover across samples * mean_PercentCover_I: Percent cover converted to an integer value * BioUnit: The bioregional unit associated with the site location * MPAorRef:  Either the name of the MPA the site was collected in or a "Ref" area associated with an MPA * depth: Depth of the site (m) * surfacereliefratio_n3: measure of surface relief using GA-SaMMT v1.2 with a neighbourhood of 3 cells (no units) * slopeposition_n10: Relative elevation of a cell compared to the overall seascape calculated using the benthic terrain modeller with a neighbourhood of 10 cells (no units) * slopeposition_n20: Relative elevation of a cell compared to the overall seascape calculated using the benthic terrain modeller with a neighbourhood of 20 cells (no units) * slopeposition_n100: Relative elevation of a cell compared to the overall seascape calculated using the benthic terrain modeller with a neighbourhood of 100 cells (no units) * surfacereliefratio_n10: measure of surface relief using GA-SaMMT v1.2 with a neighbourhood of 10 cells (no units) * surfacereliefratio_n20: measure of surface relief using GA-SaMMT v1.2 with a neighbourhood of 20 cells (no units) * surfacereliefratio_n250: measure of surface relief using GA-SaMMT v1.2 with a neighbourhood of 250 cells (no units) * surfacereliefratio_n150: measure of surface relief using GA-SaMMT v1.2 with a neighbourhood of 150 cells (no units) * surfacereliefratio_n50: measure of surface relief using GA-SaMMT v1.2 with a neighbourhood of 150 cells (no units) * surfacereliefratio_n5: measure of surface relief using GA-SaMMT v1.2 with a neighbourhood of 5 cells (no units) * slopeposition_n50: Relative elevation of a cell compared to the overall seascape calculated using the benthic terrain modeller with a neighbourhood of 50 cells (no units) * slopeposition_n250: Relative elevation of a cell compared to the overall seascape calculated using the benthic terrain modeller with a neighbourhood of 250 cells (no units) * slopeposition_n500: Relative elevation of a cell compared to the overall seascape calculated using the benthic terrain modeller with a neighbourhood of 500 cells (no units) * surfacereliefratio_n500: measure of surface relief using GA-SaMMT v1.2 with a neighbourhood of 500 cells (no units) * landform_n3: Morphology of the seafloor based on the size, slope, length-width ratio using GA-SaMMT v1.2 with a neighbourhood of 3 cells (no units) * landform_n100: Morphology of the seafloor based on the size, slope, length-width ratio using GA-SaMMT v1.2 with a neighbourhood of 100 cells (no units) * landform_n250: Morphology of the seafloor based on the size, slope, length-width ratio using GA-SaMMT v1.2 with a neighbourhood of 250 cells (no units) * landform_n500: Morphology of the seafloor based on the size, slope, length-width ratio using GA-SaMMT v1.2 with a neighbourhood of 500 cells (no units) * landform_n10: Morphology of the seafloor based on the size, slope, length-width ratio using GA-SaMMT v1.2 with a neighbourhood of 10 cells (no units) * landform_n50: Morphology of the seafloor based on the size, slope, length-width ratio using GA-SaMMT v1.2 with a neighbourhood of 50 cells (no units) * roughness_n10: A measure of seafloor complexity calculated using the variation in slope and aspect vectors calculated using Benthic Terrain Modeller with a neighbourhood of 10 cells (no units) * roughness_n250: A measure of seafloor complexity calculated using the variation in slope and aspect vectors calculated using Benthic Terrain Modeller with a neighbourhood of 250 cells (no units) * roughness_n50: A measure of seafloor complexity calculated using the variation in slope and aspect vectors calculated using Benthic Terrain Modeller with a neighbourhood of 50 cells (no units) * roughness_n100: A measure of seafloor complexity calculated using the variation in slope and aspect vectors calculated using Benthic Terrain Modeller with a neighbourhood of 100 cells (no units) * roughness_n20: A measure of seafloor complexity calculated using the variation in slope and aspect vectors calculated using Benthic Terrain Modeller with a neighbourhood of 20 cells (no units) * s_sst_c: Average summer sea surface temperature. SST values are associated with the year of survey (degrees Celsius) * a_acs: Average annual current speed. Current speed values are associated with the year of survey (m/s) * w_mov: Average winter wave orbital velocities. Orbital velocity values are associated with the year of survey (m/s) * LocalRetention: The proportion of released spores that settle back to the source site from the biophysical modelling. Values are associated with the year of survey (no units) * InDegree: The total number of significant connections coming into a site from the biophysical modelling. Values are associated with the year of survey (no units) * Easting_VG94: Easting value from the VICGRID 1994 coordinate system (m) * Northing_VG94: Northing value from the VICGRID 1994 coordinate system (m) ## Code/software The R statistical software was used for all analyses. The mgcv package was used for all GAMMs models (v 1.9-1) Predictions from the GAMMs were done using the terra package (v 1.7-46) The gratia package was used to make the partial dependence plots (v 0.8.1) ## Access information Other publicly accessible locations of the data: * bathymetry maps: [https://portal.ga.gov.au/persona/marine](https://portal.ga.gov.au/persona/marine), [https://mapshare.vic.gov.au/coastkit/](https://mapshare.vic.gov.au/coastkit/) * [https://portal.aodn.org.au/](https://portal.aodn.org.au/) * [https://zenodo.org/](https://zenodo.org/)