Description of the Trial Site for Experiments 1 and 2

The trial site was located alongside Highway 5, south of Taupo near the Rangitaiki River, New Zealand (176.34°E, 38.86°S; 962 m above sea level) and consisted of scattered wilding P. contorta in the height range of 2 to 8 m. Trees were aged at between 3 and 9 yr. The region is commonly referred to as the “frost flats” due to the common, unseasonal frosts and a cool, humid climate (mean annual rainfall of 1200 to 1600 mm, midsummer mean temperature of 14.5 C and midwinter mean temperature of 4.1 C). Besides the exotic conifers, the predominant vegetation at the site was evergreen sclerophyllous shrubland comprised mostly of the epacridaceous shrub monoao (Dracophyllum subulatum Hok. F.) with an understory of lichens and moss. The Rangitaiki region, or southeastern corner of the upper Kaingaroa Plateau, consists mainly of undulating plains overlain by Taupo Pumice colluvium 0.6- to 1.0-m deep. Underlying the Taupo Pumice colluvium is ca. 6 to 8 m of earlier rhyolitic tephras, particularly Waimihia Ash (Smale 1990).

Experiment 1

Seventy-five trees in the height range of 3 to 5 m were selected for use in the first dose–response trial. The trees were carefully selected to ensure they were at least 5 m from adjacent trees. The height (HT), mean tree crown (canopy) diameter (CD), average of the crown diameter at the point of maximum crown width and the diameter perpendicular to this line, and diameter at breast height (DBH at 1.4 m) of each tree was measured before treatment application (Table 1). Following measurement, the trees were ranked from smallest to largest and randomly allocated to treatments to ensure an even distribution of heights across treatments. A product containing 12% triclopyr in an MSO (120 g L^-1 triclopyr butoxyethyl ester in ethyl esters of fatty acids derived from canola oil; X-Tree^® Wet and Dry, United Phosphorus Limited, Auckland, New Zealand), and widely used in operational spot control of wilding conifers, was then applied using a purpose-built bark-spraying unit (Drench Gun Services, Palmerston North, New Zealand;  www.drenchgunshop.co.nz). The product was applied to the whole stem of selected trees in February 2018. The spray equipment was fit with an adjustable flat-fan nozzle and could be calibrated to deliver specific volumes using 20-ml aliquots. Four doses (Dose) of triclopyr (9.6, 19.2, 30, and 60 g) were applied as four product volumes (80, 160, 250, and 500 ml) on each of 15 trees, with an additional 15 trees included as untreated controls (Table 1). Maximum and minimum temperatures on the day of herbicide application were 23.7 and 12.9 C, respectively. Twelve months after treatment, all trees were assessed for signs of herbicide damage where the upper and lower halves of each crown were scored for percentage needle browning in 10% increments. In alignment with operational treatments for an acceptable result, we chose 90% as the threshold for an acceptable level of control and defined this as a mortality threshold, while recognizing that some trees may survive. Trees within the 90% to 100% brown foliage class were therefore categorized as “dead.”

Table 1.

Treatments applied to the whole stem of Pinus contorta trees (n = 15 per treatment) and tree size attributes in Experiment 1 with the mean and standard error for the tree size variables height (HT), crown diameter (CD), and diameter at breast height (DBH) shown.

Experiment 2

A second bark application trial was carried out at the same site in February 2019. A range of doses (Dose), including lower levels of triclopyr per meter height (with the same product) than were used in Experiment 1, were applied to the bark of trees. Application method (Method) was either over the whole stem or as per label recommendations, which is a band of product applied from the base of the tree to a stem height equivalent to three times the DBH of the trunk (Anonymous 2017). The dose was calibrated as active ingredient (ai) delivered per unit of tree height with five doses tested: 0, 1.5, 3, 4.5, or 6 g triclopyr applied m^-1 tree height (equivalent to 12.5, 25, 37.5, and 50 ml m^-1 tree height, respectively). Trees were measured before treatment allocation. Following measurement, trees were ranked according to size and then allocated to one of five treatments, such that each treatment included equivalent size range and exact dose per tree could be calculated. Each treatment was applied to 15 trees in the size range ∼2-to 8-m height (Table 2). Maximum and minimum temperatures on the day of application were 20.5 and 10.4 C, respectively. Twelve months after treatment, trees were assessed in the same manner as for Experiment 1. Similarly, trees with 90% to 100% brown foliage were classified as dead.

Using the bark-spraying unit (same as for Experiment 1), all applications were made to the bark on the main stem of the tree, starting at the base of the tree, moving upward, and either covering the whole stem or limited to the base of the tree only, as per treatment allocation assigned to the tree.

Table 2.

Treatments applied to Pinus contorta trees (n = 15 per treatment) and tree size attributes in Experiment 2 to determine the effect of application method and herbicide dose on tree mortality with the mean and standard error for the tree size variables height (HT), diameter at breast height (DBH) and crown diameter (CD) shown.

Operational Spraying and Survey Data

Operational (aerial) spot spraying using the ABBA technique was carried out by the Department of Conservation in the Kaweka Conservation Park (176.27°E, 39.46°S). The control area covered 4,000 ha and was carried out in the summer of 2016 and 2017. Triclopyr butoxy ethyl ester, applied as X-Tree^® Wet and Dry, was applied to individual trees as per best practice (Russell 2019). Data on tree location and amount of herbicide applied were captured at the time of application for each tree using an onboard bespoke turbine-driven flow meter and GPS unit (Howell and Cockburn 2016). While the accuracy of this system has not been quantified, these flow meters are generally accurate (<1% error; Smith 2001) with calibration.

In late winter 2018, approximately 18 mo after treatment, nine belt transects (10-m wide) with a total length of 2 km were randomly placed across an area that was part of the controlled operation in 2016 to 2017 (Figure 1). Wilding conifers found in the transect area were georeferenced (accuracy >20 cm postprocessed) and their HT, crown area (CA based on two perpendicular CDs), status (dead/alive), and canopy health (% brown and green foliage) were recorded. As with the experiments, trees with more than 90% brown foliage were classed as dead. Operational spray rates and ground-based tree information were merged by matching the locations of georeferenced wilding conifers on the ground with operational locations recorded by the GPS unit and flow meter during spraying.

Analyses

Because the goal of all assessments was to determine treatments that cause mortality, all models used data in which trees were assessed as 0 = alive and 1 = dead. All analyses were conducted in R (v. 3.6.1; R Development Core Team 2014). Optimal binomial models of mortality were developed using backward selection applied to a generalized linear model (GLM), assuming a binomial distribution and logit link. Where complete separation errors occurred, the binomial GLM was fit by penalizing the likelihood using the method proposed by Firth (1993): bias reduction of maximum-likelihood estimates. All models took the general form:

where p is the probability of mortality, β₀ is an intercept, β₁‚ . . . ‚ β_p are a set of regression coefficients, and x₁‚ . . . ‚ x_p are the set of predictors and their interaction terms, k = 1‚ . . . ‚ p, typically estimated by the maximum-likelihood method (Hosmer et al. 2013).

For all experiments and the survey, trees were scored as dead when foliage browning exceeded 90%. Initial model selection was based on the Akaike's information criterion (AIC) (Akaike 1998), with goodness of fit evaluated using the likelihood ratio test (Hosmer et al. 2013) and Hosmer-Lemeshow test (H&L) (Hosmer and Lemeshow 1980). Mcfadden's generalized R² (Mcfadden 1974) and k-fold cross validation (James et al. 2013) were used as measures of predictive power. In k-fold cross validation, the original sample is randomly partitioned into k equal-sized subsamples. Of the k subsamples, a single subsample is retained as the validation data for testing the model, and the remaining k - 1 subsamples are used as training data. The cross-validation process is then repeated k times, with each of the k subsamples used exactly once as the validation data. The k results can then be averaged to produce a single estimation. The advantage of this method over repeated random subsampling is that all observations are used for both training and validation, and each observation is used for validation exactly once.

For Experiment 1 (Table 1), the full model (Equation 2) included herbicide dose (Dose; m^-1 of tree HT) and tree size (HT, DBH, or CD) as explanatory variables:

For Experiment 2 (Table 2), the full model included both herbicide dose (Dose; m^-1 of tree HT) and application method (Method: Whole Stem or Label), and their interactions (Equation 3). During model selection, the effects of the tree size variates HT, DBH, and CD on mortality across treatments were also evaluated. The significance of factor terms was followed up by applying a multiple-comparison procedure using Tukey contrast (R package EMMEANS; Lenth 2020).

The data from Experiment 1 and Experiment 2 (Whole Stem data only) were then combined and used to evaluate the best predictors of dose–response by expressing dose as a function of unit tree size, calculated as grams active ingredient of a range of tree size measures as follows:

Figure 1.

Location of the belt transects used to locate treated trees in the Kaweka Conservation Park. Red points indicate dead and green points indicate live trees.

Six models were fit by penalizing the likelihood using the bias reduction method proposed by Firth (1993). All models were evaluated using log likelihood, AIC, AICc (the corrected AIC which can be used as an approximation of AIC for a correction for small sample sizes), and the Bayesian information criteria (BIC) as criteria of goodness of fit, with predictive accuracy evaluated using a k-fold cross-validation approach.

Operational Spraying in Kaweka Conservation Lands

The relationship between the total amount of triclopyr applied to trees assessed as dead and tree size variates (HT and CD) was summarized using scatter plots and trend lines. These data were compared with data from Experiments 1 and 2 to indicate whether reductions in herbicide use could be achieved during operational programs in which triclopyr in oil is used to kill trees using the ABBA method.

Table 3.

Summary of the binomial generalized linear model analysis of mortality of Pinus contorta trees as a function of triclopyr dose (Dose), method of application (Method), and tree size (CD, HT, and DBH) (cf. Equation 3) with the regression coefficients (β) and their standard error (SE), the significance of the terms in the model (deviance, df, P-value), goodness-of-fit statistics (LRT, H&L), and statistics measuring predictive power (McFadden R², AIC) shown.^a

Figure 2.

Probability of mortality of Pinus contorta trees in the 9.6 g of triclopyr treatment (n = 15; Experiment 1). Mortality is shown both as a function of (A) tree height and (B) as grams of triclopyr applied per meter of tree height. Trend lines are indicative only, and no statistical significance is attached to these lines

Experiments 1 and 2

Barring the untreated control, for which there was no mortality, 93% of the treated trees died in Experiment 1. All surviving trees (n = 3) occurred in the lowest dose treatment in which 9.6 g triclopyr was applied in 80 ml of product (Table 1). The high mortality precluded fitting a model that included tree size and dose parameters using binomial GLMs (cf. Equation 2), as there was complete separation within the predictor variable Dose. This outcome reflected the separation in mortality between the treated (92% dead) and untreated trees (100% alive). To overcome the separation error in the predictor variable Dose, an alternative method proposed by Firth (1993) was used in the model. However, the resulting model was not significant (α > 0.05; data and model not shown) and the analysis of the data from Experiment 1 did not provide conclusive results for the effect of tree size on dose response for bark applications of triclopyr in oil for the range of concentrations and tree sizes tested. The results indicated that the lowest dose of 9.6 g triclopyr, or 2 g triclopyr m^-1 of tree height, was potentially sufficient to kill trees smaller than 4 m in height, but this dose was insufficient for larger trees (Figure 2). These results were used to inform a follow-up experiment that included a different range of dose treatments and a wider range of tree sizes.

In Experiment 2, the independent variables Dose, Method, and tree CD, and their interactions, were shown to have a significant (α = 0.05) effect on tree mortality (Table 3; Figure 3). All three tree size covariates (HT, DBH, and CD) were tested as predictors of mortality; however, the model with CD as a covariate to Dose and Method had the lowest AIC value and was therefore selected as the best model (Table 3). Herbicide application to the whole stem resulted in significantly (α = 0.05) higher mortality (90 ± 4%) than where herbicide was applied as recommended by the label (54 ± 4%), where no mortality was predicted across all tree sizes and doses tested (Tukey contrast; P < 0.05; Figure 3).

Figure 3.

Probability of mortality for Pinus contorta trees with crown diameters (CD) 1.5, 2.5, 3.5, and 4.5 m treated with 0, 1.5, 3.0, 4.5, and 6.0 g of triclopyr (applied as 12% active ingredient in methylated seed oil) applied per meter of tree height (HT) either to (A) the bark of the whole main stem of the tree or (B) as per the label recommendation (see Table 4). Equations for whole stem: logit(p_i) = -1.7601 - 0.008486*CD þ 4.64759*Dose - 0.00788*CD*Dose; and the label method: logit(p_i) = -0.8177 - 0.006517*CD þ 0.26995*Dose þ 0.00077*CD*Dose.

The lower efficacy when treating the base of the tree alone was a surprising outcome, as this method of application has long been an accepted method for control of unwanted invasive woody species (DiTomaso and Kyser 2007; Nelson et al. 2006), including for control of individual, invasive pines in New Zealand (Russell 2019). While these results may not necessarily reflect the efficacy of the aerial bark applications (ABBA method), in which most of the stem is treated with oil, they could have important implications for ground-based operations. Basal bark treatments have traditionally been most reliable on smaller stems and are known to be less effective on larger stems with thicker bark layers (DiTomaso and Kyser 2007; Schneider 1991; Yeiser et al. 1989). Variation in the efficacy of basal bark applications on wilding P. contorta has also been reported by operational staff in New Zealand, but there are no quantitative data on tree mortality from routine ground-based basal bark operations. The implications for control of wilding P. contorta, currently the most widespread and established wilding conifer in New Zealand, are either to develop methods to increase uptake and translocation of active ingredient when applied to the bark (e.g., changing herbicide concentrations or additives) or to treat a greater area of the whole stem where possible.

Prediction of Dose Response

One of the objectives of this research was to characterize the tree size and dose–response relationship for bark applications of triclopyr. The results from Experiment 2 indicated CD as a significant (α = 0.05) covariate in the prediction of efficacy when dose was applied per unit of tree height, with tree HT and DBH shown to be poorer covariates (Table 3). To complete the analysis, data from Experiment 1 and Experiment 2 (whole-stem treatment data only) were combined (n = 144) and used to evaluate the best predictors of herbicide dose. To simplify models, dose was evaluated per unit tree size expressed as HT (with and without CD as covariate), DBH, CD, CA, and CVI (Table 4;  Supplementary Table S1). Out of the six models fit, Model ii, where dose was expressed as a function of HT only, provided the poorest fit to the data (highest AIC) with the lowest predictive accuracy (lowest CV). Model vi, where dose was expressed as a function of unit tree CVI, provided the best fit, with Models i, iii, iv, and v also showing reasonable fit and predictive accuracy.

Table 4.

Model comparison for the binomial generalized linear models (GLMs) of the six models obtained by penalizing the likelihood using the bias reduction method proposed by Firth (1993), with Akaike information criterion (AIC) selected model shown in bold.^a

Table 5.

The estimated rate of triclopyr (g) per unit of tree size to incur 95% or 99% probability of mortality and associated confidence intervals (CI) for models i – v.

Although models i, iii, iv, v, and vi were not nested, we used the informal deviance tests as a rough guide for comparing models (Yee 2010). Models iii, iv, v, and vi, which were simpler models, were subsequently compared with Model i, which was a more complex model (combining amount of triclopyr per unit tree height with CD as a covariate). The approximate chi-square test statistic under the null hypothesis that Model iii was comparable to Model i yielded a P-value of Pr(χ²₂ >19:92) < 0:001, indicating that Model i was a better fit compared with Model iii. The approximate chi-square test statistic under the null hypothesis that Model iv was comparable to Model i yielded a P-value of Pr(χ²₂ >5:15) = 0:076, indicating that Model iv was a better fit than model i. The approximate chi-square test statistic under the null hypothesis that Model v was comparable to Model i yielded a P-value of Pr(χ²₂ >9:59) < 0:05, indicating that Model i was a better fit compared to Model v. The deviance for Model vi was 48.1, while that for Model i was 45.3. A comparison yielded a χ² of 2.77 on 2 degrees of freedom, suggesting that Model vi, too, appeared a better fit compared with Model i. However, based on goodness of fit and prediction accuracy models i, iii, iv, v, and vi all provided a good fit to the data with reasonable predictive accuracy, with model vi > iv > I > v > iii > ii.

Using models i to vi, the dose of triclopyr to be applied to the whole stem and estimated to cause mortality with 0.95 or 0.99 probability was predicted (Table 5;  Supplementary Figures S1–S5) and used as a basis to determine the total amount of active ingredient required to kill 10 randomly selected trees from the trial data set (Table 6). The purpose of this analysis was to evaluate the sensitivity of model choice to the estimate of the amount of herbicide required to kill a randomly selected set of trees of differing sizes. These models varied in prediction of effective dose by more than 50% for any one tree. For example, Model vi predicted that 6.2 g of triclopyr (or 52 ml of product) would be required to kill a tree (probability = 99%) with a HT of 3.3 m and a CD of 1.9 m, while Model ii predicted 17.0 g of triclopyr (or 141 ml of product) for the same tree. However, while these calculations indicated that the six models predicted different amounts of active ingredient to kill individual trees, using either the 95% or 99% probability values, the total amount of product required to kill all 10 trees was similar. Applying these models over the whole data set gave similar predictions on the amount of product required. As these trials were conducted at only one site in New Zealand, ideally, validation trials would be conducted to determine the application of these models to a wider array of sites and species.

Other studies have shown similar effects of tree size on treatment efficacy in relation to dose (Enloe et al. 2016). Crown diameter, which is strongly linked to crown area and crown volume, has been shown to be a good predictor of herbicide treatment efficacy (Knowe 1991; Miranda-Fuentes et al. 2016), no doubt because of the strong relationship of this variate to woody plant biomass (McGinnis et al. 2010). In a study to develop equations to predict woody plant biomass using nondestructive assessments, McGinnis et al. (2010) found CD, crown volume, and both CD and height as the best predictors of biomass (albeit for broadleaf woody species). In a pine release study, Zedaker and Freyman (1988) reported the highest control of hardwoods in relation to herbicide applied was based on calculated crown volume, intermediate control from crown area, and the lowest control calculated from crown height. Similar outcomes were also reported by Knowe (1991). Expressing dose per unit of DBH was also shown to be a good predictor of efficacy in the combined data set, an unsurprising result given the strong relationship of this variable to tree biomass (Ung et al. 2008). In operational practice, however, CA, CVI, and DBH may not be practical variables to estimate visually and are possibly unsuitable for dose estimation ad hoc (e.g., McGinnis et al. 2010). While DBH may be an easy variable to estimate in forest stands, for individual or clumped weedy trees, the stem is often obscured by foliage, with multiple stems not uncommon. In future, application of herbicides from a UAV to either the crown foliage or bark might allow for estimation of crown volume or crown area in real time, allowing for more precise application of spray doses calibrated to two-dimensional or three-dimensional estimates of tree size. In the meantime, operators applying herbicide to the tree from the ground or a helicopter are more likely to use size variates easier to assess visually, such as tree HT or CD. Out of the size variates HT and CD, the models developed in this study indicated that dose applied on an estimate of CD provided a reasonable fit with good prediction accuracy.

Table 6.

Total triclopyr (g) in oil predicted to provide a 95% and 99% probability of mortality of 10 randomly selected trees when applied to the bark.^a

Figure 4.

Data collected during an operational aerial spot-spraying exercise using the aerial basal bark application, or ABBA, method of applying triclopyr in methylated seed oil to individual conifers. Shown are the amounts of herbicide applied to trees that were evaluated as dead (x) or alive (o) expressed as the total amount of triclopyr (g active ingredient) applied to trees of different (A) height and (B) crown diameters. Gray lines show the trend (best fit) of the operational data. Shown as dotted red lines on each figure are also the comparative rates of triclopyr shown to kill trees in a research trial where herbicide was applied direct to the bark of the trees.

Operational Spraying and Survey

A total of 460 trees were assessed and spatially georeferenced across the nine sampling transects, with 270 of these aligning to referenced points in the operational data set. Of the remaining 190 trees, 47 could not be assigned to a recorded location and herbicide amount, and 143 trees were recorded as not being sprayed, that is, 31% of trees assessed were missed during the spraying operation. These untreated trees ranged in HT from 0.3 to 5.3 m and followed a similar distribution in HT to those trees that were treated. The high percentage of trees missed is most likely a reflection of the steep terrain and high density of surrounding arboreal shrubs such as mānuka [Leptospermum scoparium (J.R. Forst. & G. Forst.)] and rawirinui [Kunzea robusta (de Lange & Toelken)], which would have reduced the visibility of conifers during the aerial spraying operation.

The 270 trees aligning to treated points ranged in size from 0.5 to 6 m, with an average HT and CD (±SE) of 1.9 ± 0.07 and 1.0 ± 0.03 m, respectively. Of the trees treated, 86% (231) were killed. Trees that were not killed ranged in height between 0.5 and 3.9 m.

The relationship between tree size (HT and CD) and the amount of herbicide applied was plotted for all trees aligning to treated points (Figure 4). The scatter in the operational data highlighted the wide range in the amount of herbicide applied to trees of similar size. When tree size was described using either HT or CD, the lower bounds of the operational rate range that achieved mortality were very similar to the predicted rates derived from Experiments 1 and 2 (Figure 4). This result implies that in many cases, operational herbicide rates can be up to about seven to eight times higher than rates required to achieve mortality. However, it is notable that the doses applied to trees not killed were mostly between the modeled dose and the center of the distribution, highlighting that such threshold doses may not be optimal for application from a helicopter, probably due to application errors and a significant proportion of spray missing the target tree.

Herbicide application methods that involve treating individual stems with herbicides can result in large amounts of herbicide being used if applied to scattered individuals. Overdosing occurs when practitioners exceed herbicide label rates or do not follow area application rates, an activity that is easy to do in conservation lands that are not delineated or measured as for agricultural or forestry land use (Holmes and Berry 2009). Of all the actions to reduce the environmental impact of herbicides, dose adjustment is one of the most important. As a first step, the data generated from the experimental trials could be used to provide better recommendations to ground operators on the amounts of product to apply to trees of different sizes. This would reduce both the cost of the ground operations and input of herbicide into the environment. The data from the aerial control operation carried out in the Kaweka Conservation Park indicated a high variability in the amount of herbicide applied to trees of similar sizes. Considering the difficulties of targeting and tree size estimation from a hovering helicopter, this is not at all surprising. From an operational perspective, the additional costs of excess herbicide are probably negligible considering the cost of flying, personnel time, and risk of managing to a minimum dosage threshold that may require repeat operations. However, a comparison of the amount of herbicide used operationally with that estimated to provide 95% probability of mortality from experimental trials indicated that dose reductions by at least 50% were feasible. Hence, rate reductions and improved targeting of dose to tree size could be encouraged from an environmental perspective, even if economic gains are marginal.

It is interesting to put a value (31%) to the percentage of trees that were missed in the ABBA operation carried out in the Kaweka Conservation Park. This is the first assessment to benchmark the efficiency of the ABBA method in steep, hilly terrain where wildings are emerging within regenerating forest. The Kaweka range is an environment with flat hilltops dominated by subalpine shrub-lands and tussocks, steep hill flanks covered in rapidly regenerating L. scoparium and kānuka [Kunzea ericoides (A. Rich.) J. Thomp.], and gullies comprising shrub species mixed with beech (Nothofagaceae Kuprian) and indigenous conifer species (Podocarpaceae Endl.). Not only is this difficult terrain to navigate by helicopter, but this is also a difficult vegetation type from which to identify the emerging wilding conifers. A value of 31% missed is reasonably high and highlights a requirement for a repeat operation with the associated high helicopter and personnel costs. Recent advances in detecting wildings remotely (Dash et al. 2019b) confirmed that in simpler landscapes (grasslands), trees taller than 1 m can be detected with high accuracy, and with further improvements in the technology (e.g., UAVs combined with multispectral sensors) and algorithms (Dash et al. 2019a), it is expected that even in complex landscapes, the detection levels can be improved to avoid frequent control operations. Automated mapping and estimation of tree sizes based on high-resolution imagery can also provide valuable information to inform operational managers of where the greatest gains in efficiency of control operations could be made, in conjunction with better targeting of herbicide dose to tree size.

This study has provided a first step in quantifying the triclopyr dose required to kill individual P. contorta trees when applied to the bark. While the experimental work was repeated over 2 yr, the next step will be to validate these initial models across a broader range of sites and species and to test whether season of application influences efficacy. There are currently conflicting opinions from operational managers on timing effects, but there are currently no data on this topic for quantitative analysis. Further testing will allow for refinements in models and potential reductions in the amounts of herbicide used.