A historical dataset of biofouling and ballast water declarations was compiled for vessels that visited New Zealand during 1998 to 2008, which contained records for 15,745 voyages by 2,607 vessels. This dataset contained information on the vessel’s physical and voyage characteristics, and ballast discharge volumes detailed for each ballast tank. This dataset was examined to determine: (i) unique voyage and port visits by vessels and volumes of ballast discharged per vessel at each port, and (ii) management of the ballast and intention to discharge in New Zealand. The response variables used for ballast water modelling were the port of discharge and discharge volume per port. Vessel type, ballast capacity or Deadweight Tonnage (DWT), arrival port, and the intention of discharge were candidate predictor variables for estimating how much ballast water was discharged at each port. Arrival port is the first New Zealand port a vessel encounters in a voyage. Many voyages declared more than one port of ballast discharge, i.e., each journey has a first arrival port, but the vessel might visit other ports as part of the same voyage and discharge ballast. For example, a vessel’s voyage to New Zealand may include visits to Whangarei, Tauranga and Auckland, with intended ballast discharge at Whangarei and Tauranga but not Auckland; here the recorded intended ballast discharge volumes for the first two ports (i.e., those at which ballast was discharged) were tallied. The historical data contained 15,628 unique voyages out of which 1,196 journeys had more than one port of ballast discharge. Out of 17,703 stops at ports of arrival, there were 6,287 incidents of reported ballast discharge. Note that the discharge may have presented an acceptable biosecurity risk, if the discharge were carried out in a regulated way, for example after open-ocean exchange. ‘Intent to discharge’ was not used as a predictor because there was poor agreement between the declared intention to discharge and the reported discharge volume. More information about this agreement can be found in the Suppl. material 1: Table S1 in the SM summarises the description of the variables used in the historical dataset for the ballast water analysis.

Biofouling on the submerged surfaces of 322 international merchant vessels was sampled by Inglis et al. (2010) upon arrival to New Zealand. This included 166 container and general cargo vessels, 49 passenger vessels, 37 roll-on roll-off (RoRo) / car carriers, 31 bulk carriers, 21refrigerated cargo ships, 12 tankers and 6 vessels of other miscellaneous class types. A standardised sampling design was used to collect samples of biofouling from the hull and external niche area surfaces of the vessels. A questionnaire administered to the vessel’s master provided information on the design features of the vessel, hull maintenance record (e.g., age and type of antifouling coatings), and recent voyage history.

Data on the presence and mass of biofouling from 314 of the merchant vessels and associated data from the questionnaire were used to build the biofouling model. As the contemporary arrival data to which the model would predict contained only limited information on the recent voyage history of the vessels, the models were fitted using variables from the historic study that could be easily matched to contemporary information collected by BNZ. They included design specifications of the vessel (e.g., class type, design speed), vessel age, niche area, the age of antifouling coatings on arrival, and the maximum periods of inactivity (lay-up) in any location. For biofouling data, the main two response variables were biomass of biofouling per square metre and presence/absence of any biofouling. The latter had a value of 1 whenever biofouling was present, and zero if there were no biofouling. This binary variable was used to correct mass; that is, whenever mass had a zero value, but biofouling was present, a correction value of 0.001 g was added to the mass. Predictor variables in the historical biofouling data were the age of antifouling coating, vessel type, vessel age, vessel speed, season, niche, and maximum duration of lay-up. The description of the variables and their ranges are given in the Suppl. material 1: Table S2.

Contemporary data sourced from BNZ on international vessel arrivals into New Zealand were a combination of two data sets. The first included data on vessel arrivals between January 2015 to December 2017, and was derived from mandatory pre-arrival documentation provided by the vessel master to BNZ (Ministry for Primary Industries 2019). These data were originally collected to aid in the tactical evaluation of biosecurity risks. The second data set included port-by-port arrivals by merchant vessels, including the ports of first arrival to New Zealand and coastwise voyages for the year 2016. These data were purchased from the Lloyd’s Maritime Intelligence Unit, and used to verify the completeness of the BNZ data. Details of all data preparation for the analysis are provided in the Suppl. material 1 – part 1. The same contemporary data were used for both ballast water and biofouling analysis, but a different set of predictors were selected for each.

The ballast discharge model was used to predict the total annual volume of ballast water discharged at each New Zealand port. The contemporary data lacked ballast discharge information, so we had to predict it from historical data. The models were built to predict ballast discharge at each port as a function of arrival port and ship characteristics.

The historical biofouling data from Inglis et al. (2010) were used to fit a model that predicted the presence and biomass of biofouling on vessels that arrived in New Zealand between 2015 and 2017. The total mass of biofouling on a vessel was estimated as the product of its estimated mass per unit niche area (i.e., density) and the total niche area for the vessel. A model was constructed from historical data to predict the presence and per unit area of biofouling in niche areas of the vessels as a function of vessel characteristics such as the age of antifouling coatings on the vessel, and the frequency and duration of lay-up periods. As most biofouling occurs in niche areas (Inglis et al. 2010; Davidson et al. 2016), the total mass of biofouling on a vessel was estimated by calculating the total area of niches on each vessel using relationships developed by Moser et al. (2016) and scaling the predicted density for the entire vessel hull area. This model was then applied to the contemporary vessel data to predict port-level exposure to biofouling. Additional details about the modelling approach can be found in the Suppl. material 1 – part 3.

This section reports construction of an algorithm to predict the volume (and at which ports) a vessel will discharge their ballast water, given their vessel-level characteristics (e.g., type, DWT, Ballast capacity) and journey-level information. The complex aspect of this endeavour is that after entering New Zealand’s waters, each vessel may visit more than one port, with the possibility of discharging ballast water at one or more of the ports visited. As a result of discharging at one or more ports, the variable of interest is structured as a vector for each vessel, so any algorithms investigated are required to take this vector-based outcome into account.

An algorithm based on k-Nearest Neighbours (KNN, Fix and Hodges 1951) was developed to simultaneously predict which ports a vessel may visit and the subsequent discharge amount. For each journey in the contemporary data, this algorithm finds the k most similar journeys in the historical data based on some measure of similarity and uses the average ballast discharge from those journeys as a prediction of the discharge for the contemporary journey.

Denoted by X the n × p matrix of data used for model training (the training data is the historical data as described in the Suppl. material 1 – part 1), where n is the number of voyages (observations) and p is the number of fields (predictors) in the data; let Y be the corresponding n × v matrix of discharge volumes, where v is the number of ports. Let z_j be a vector of length p, containing the predictors for the j^th voyage in the contemporary data, for which we wish to predict ballast discharge,ŷ_j, at each of the ports. KNN algorithm as applied to ballast discharge data can be explained as below:

1. Repeat for j = 1, …, m, where m is the number of voyages in the contemporary data for which we wish to predict:

1. Calculate D (z _j , X), the distance from z _j to each of the rows xi, i = 1, …, n in the training matrix X.
2. Find S, the set of the k ^th smallest distance in D (if there are ties, keep these in the set).
3. Calculate

${\hat{y}}_j=\frac{1}{\left|S\right|}\sum _{i\in S}{y}_i$

i.e., the predicted discharge values are the average of the k^th nearest neighbours discharges.

1. For each year t in the contemporary data:

1. Find S _t[v], the set of all voyages arriving at port v during year t.
2. Calculate

${\hat{B}}_{t\left[v\right]}=\sum _{i\in S_{t\left[v\right]}}{\hat{y}}_{t\left[j\right]}$

i.e., B̂_t[v] is the total predicted ballast discharge in port v during year t.

The KNN algorithm does not work for categorical predictors, as there is no defined distance metric for mixtures of categorical and numerical predictors. To allow for categorical predictors (such as vessel type), the KNN algorithm was implemented on numerical predictor matrices created within each level of a category. For example, to use vessel type and dead weight tonnage as predictors, the algorithm was run on the matrix X formed from the subset of voyages from each level of vessel type until all predictions were made.

The historical data were structured such that each row of the dataset represented a single event in a vessel’s voyage. For example, one row may consist of the arrival event, with no discharge, and the next row a discharge event for the same vessel at a different port. It was important that if a vessel did not discharge in a port, then this non-event was not recorded. In order to have a complete vector response for each voyage, it was reasonably assumed that non-events resulted in a 0 discharge being recorded. An example of this transformation is provided in Suppl. material 1: Fig. S2 in the Suppl. material 1 – part 4.

The choice of which predictors to include, and how many nearest neighbours to use in the prediction was made using cross-validation. Given our prediction target of total ballast water volume at a port within a year, our cross-validation strategy was to leave out a single year at a time to form the testing set, with the remaining year’s data forming the training set. A schematic of cross-validation used for KNN analysis is given in Suppl. material 1: Fig. S3 in the Suppl. material 1 – part 4.

In order to provide an unbiased validation of the chosen model and help understand its operational characteristics, the data were split into two: model-fitting data comprising all voyages within the years 1999 to 2005 (inclusive) and validation data comprising all voyages within the years 2006 and 2007.

To estimate the prediction error from the KNN algorithm, a bootstrap procedure was used. Voyages were randomly sampled with replacement from the training data, and the KNN algorithm fit on each bootstrap sample. The standard deviation (and mean/median) of total discharge within a port was calculated from 100 bootstrap resampled values of total discharge.

Several different modeling approaches were applied to the historical data in the process of model selection to estimate biofouling biomass for the contemporary dataset (Suppl. material 1 – part 6). Random forests provided the best model to obtain the prediction (Breiman 2001; Liaw and Wiener 2002; see Hastie et al. 2017 for a discussion) and we used quantile-regression random forests to obtain the distribution to represent the prediction uncertainty (Meinshausen 2017).

Random forests analysis was applied to the data because the algorithm is indifferent to the conditional distribution of the response variable, meaning that no special provisions need to be taken to allow for the heavy zero inflation. The random forest model explained about 15% of the variation in the biofouling mass, which was better than its alternatives. Quantile random forests was used to obtain journey-specific prediction distributions which could then be aggregated to obtain overall levels of uncertainty for the exposure of each port to biofouling for each year.

The likelihood of entry via biofouling was assumed to be proportional to the biofouling exposure at each port, each year. For a single vessel, this was assumed to be the predicted total mass of biofouling in the niche areas of the vessel (not including the hull as a niche area). However, the contemporary data did not contain biofouling mass measurements for each vessel but did contain measurements for each vessel’s niche area in metres squared. To handle this, two models were constructed: a model to predict the biofouling mass per vessel as a function of several variables, and another model to predict the quantiles of the prediction distribution. The random forest model used the following continuous predictor variables: dry weight tonnage, days since last antifouling paint, niche area, and mean speed. Vessel type was assessed and was excluded from the models because the improvement of fit was insufficient.

To predict the density of biofouling for each vessel in the contemporary data, missing data were first handled using multivariate imputation via chained equations (Buuren and Groothuis-Oudshoorn 2010) to make five complete variant datasets, which were analysed in parallel. The random forest biofouling presence model was applied to each contemporary vessel movement, which gives an estimate of the quantum of biofouling. Then the prediction distributions were estimated using the quantile random forest. Once the density of biomass present on each vessel was predicted, it was multiplied by the niche area of the vessel to give a total niche biomass prediction, which was then used in the allocation procedure.

We calculated an entry likelihood score from ballast discharge and biofouling propagule pressure on a scale of 0 (lowest) to 1 (highest). That is, the mean values from the estimated prediction distribution were scaled to each of the respective discharge and biofouling amounts by the maximum amount within the ports, which sets the scaled score to 1 for the port with the highest amount. We did not subtract the smallest so that no port was allocated a weight of 0. These scores were then averaged to give the relative entry likelihood score for each port. Weights were then calculated in proportion to these relative scores. The weights sum to 1 across the 15 ports and could, therefore, be used to allocate total surveillance effort among ports relative to the likelihood of entry of marine NIS for each port.

We fitted the KNN model using all combinations of vessel type and arrival port as categorical variables, DWT or date of arrival as the nearest neighbour matrix, and k = 1,3,5 nearest neighbours to the model data set. Results of the cross-validation in Table 5 (Suppl. material 1 – part 7) showed that the KNN that grouped by both port of arrival and vessel type included DWT for calculating nearest neighbours, and using five nearest neighbours performed the best, i.e., this combination of distance variables had the smallest RMSE value (88732.8 for k = 5). A dataset of years 2006–2007 was used to test KNN validation by comparing the mean of predicted ballast water values (as predicted by the bootstrap) with the true values (Suppl. material 1: Fig. S18 in Suppl. material 1 – part 8). Based on the results of the KNN validation summarised in this Fig., the KNN models performed reasonably well for prediction; the observations follow the 1:1 line and double standard deviation lines for the port-level discharge mostly cross the line.

The best model arising from the cross-validation fitting study of the KNN algorithm as selected was used to predict contemporary ballast discharge per port. For this application, the full historical data (1999–2007) were used for training. Fig. 1 shows the mean and 10% and 90% quantiles as estimated by bootstrapping. Tauranga, with an average of 2255.73 ± 292.86 kilotons, had the highest predicted average ballast discharge across the years 2015–2017, followed by New Plymouth (1025.22 ± 231.89), Lyttelton (907.44 ± 129.39), Taharoa (831.84 ± 242.45), Napier (812.54 ± 191.47), Whangarei (624.92 ± 165.55), and Gisborne (575.33 ± 160.74). For a port receiving a high number of vessels, Auckland had a low volume of discharge by comparison. Following an increase in 2016, the predicted ballast discharge decreased in 2017 in almost all ports (Fig. 1). The lowest predicted values for ballast discharge were observed in Milford Sound, Picton and Bluff, and zero for Westport.

Figure 1.

Predicted ballast water discharge (kilotonnes) per port for the contemporary data. Data are the mean and standard deviation of total volume as estimated by bootstrapping. The figure displays the mean and 10% and 90% quantiles of the bootstrap distribution.

Based on the Random Forest statistics, e.g., Gini impurity index and MSE (mean squared error), important variables for biofouling prediction were days since last antifouling paint, niche area, vessel speed and DWT. Variability in the modelling was accounted for by predicting quantiles of the vessel biomass (Fig. 2).

Figure 2.

Predicted vessel biofouling exposure (tonnes) to each arrival port for the contemporary data (2015, 2016, 2017 and combined years). The median (dot) and 10% and 90% quantiles (upper and lower lines) of the prediction distribution are displayed.

As summarised in Table 1, Auckland had the highest average predicted biofouling exposure, followed by Tauranga, Lyttelton, Napier, and Wellington. Ports with very low biofouling exposure were Whangarei, Gisborne, Picton, and Taharoa. Fig. 3 shows the variation of the actual biofouling and the prediction distribution, which demonstrates that the model is not of particularly good quality. That is, some of the data points, particularly for higher values of biofouling mass, are well below the 1:1 line, indicating an underestimation of higher values, which reflects regression to the mean.

Figure 3.

The variation of the actual biofouling mass (g/m2) versus the predicted biofouling mass (g/m2). The vertical grey lines represent the 10% and 90% quantiles of predicted values. The diagonal line shows the 1:1 relationship.

Finally, for each arrival port, the biofouling exposure can be compared with the ballast discharge (Fig. 4). This shows that the arrival ports of Auckland and Tauranga have the highest exposure for biofouling, and Tauranga has the highest volume of ballast discharge (i.e., highest entry likelihoods). The next cluster of ports comprises New Plymouth, which had the second highest ballast discharge but received much lower biofouling exposure than many other ports; Wellington, which had low volume of ballast discharge and high exposure to biofouling; and Lyttelton and Napier, which seemed to have relatively high levels of exposure to both ballast discharge and biofouling. The balance of ports have relatively low exposure.

Figure 4.

Predicted vessel ballast water discharge (kilotonnes) vs biofouling exposure (tonnes) per port. The figure displays the mean and 10% and 90% quantiles of the bootstrap distribution for ballast water and from the prediction distribution for biofouling. The values are the average of ballast water and biofouling mass over years 2015 – 2017.

Based on the relative likelihood of entry of marine NIS via these two pathways, the highest surveillance weights would be allocated to the ports Tauranga and Auckland, followed by Lyttelton, Napier, New Plymouth, and Wellington, whereas the lowest surveillance effort would be allocated to the ports Gisborne, Timaru, Bluff, and Picton.

The relative entry likelihood weighting, calculated from the mean values of the estimated prediction distributions and scaled for the respective ballast discharge and biofouling exposure, is provided in Table 2. Based on these results, the ports with higher likelihood of marine NIS entry were Tauranga, Auckland, Lyttelton, Napier, and New Plymouth with scores of 1, 0.59, 0.43, 0.38, and 0.31, respectively. Tauranga was assigned a relative allocation weight of 0.23, followed by Auckland (0.13), Lyttelton (0.1), and Napier (0.09). New Plymouth, Wellington, Nelson, Dunedin, Whangarei, and Taharoa were ranked next based on their surveillance allocation weights. Gisborne, Timaru, Bluff, and Picton had the lowest relative weighted scores. There was a positive relationship between the number of vessel arrivals per port and the entry likelihood scores (r² = 0.75). Ports where the entry likelihood is greater than might be predicted by numbers of arrivals alone include Tauranga, Napier, New Plymouth and Taharoa, where there are relatively large proportions of bulk carriers or tankers. Ports where entry likelihood seems lower include Auckland and Wellington. To investigate the variability of the ranking of sites, this allocation was also performed to each of the bootstrap replicates. Fig. 5 shows the histogram of ranks from this procedure.

Figure 5.

Rank histograms for surveillance allocation, as estimated by bootstrap and prediction simulation.