Trade-offs between environmental and economic factors in conversion from exotic pine production to natural regeneration on erosion prone land

Background: Some of New Zealand’s exotic pine (Pinus radiata D.Don) forests were planted for erosion mitigation but cultural, legislative, environmental, and profitability limitations in some parts of the landscape have led to reassessment of their suitability. There is limited information to support landowner decisions on the viability of natural regeneration of native forest post-pine-harvest. Methods: We evaluated scenarios of post-harvest natural regeneration, compared to remaining in pine production, using erosion susceptibility determined from historical occurrence of landslides, gullies and earthflows, biophysical growth modelling of mānuka–kānuka (Leptospermum scoparium-Kunzea ericoides (A.Rich) Joy Thomps.) shrubland using the process-based CenW model, and cost-benefit analyses using NZFARM with two land use change scenarios, at two levels of erosion mitigation ± honey profits. Results: In our study area, the Gisborne Region (North Island of New Zealand), ~27% of the land has moderate–very high susceptibility to landslides, 14–22% a high probability of contributing material to waterways, and 19% moderate–very high gully erosion susceptibility. Pines grow 10 times faster than naturally regenerating mānuka–kānuka shrubland, but mānuka–kānuka is used for honey not wood production. Natural regeneration resulted in losses of $150–250 ha−1 yr−1 compared to the current profitability of pine production. Honey production offset some reduction in pine revenue, but not fully. Thus, the viability of shifting from pines to native forest is highly dependent on landowner impetus and value for non-market ecosystem services (such as cultural and biodiversity values) provided by native forest. Conclusions: A mosaic of land uses within a property may sufficiently offset income losses with other benefits, whereby highly erosion-prone land is shifted from rotational pine forest production to permanent native forest cover with honey production where possible. At the regional scale in Gisborne, the conversion of the most highly susceptible land under production forestry (315–556 ha) to natural regeneration has the potential for wider benefits for soil conservation reducing erosion by 1–2.5 t yr−1 of sediment facilitating achievement of cleaner water aspirations and habitat provision.


Introduction
Some owners of exotic pine (Pinus radiata D.Don) forests in New Zealand aspire to convert parts of their holdings to native forest using natural regeneration after pine harvest. The viability of converting from profitable pine production to natural regeneration is unclear and the lack of information on costs and opportunities associated with adopting this conversion strategy inhibits informed decision-making by landowners.
Much of New Zealand was cleared of native forest during Polynesian and European colonisation (McGlone 1983;Guild & Dudfield 2010); however, some cleared areas were highly susceptible to erosion (McGlone 1989;Rhodes 2001;Guild & Dudfield 2010). Nonnative conifers were planted by the NZ Government in the 1960s as a fast-growing solution to mitigate erosion. Government-owned forests in the Gisborne Region were privatised in the 1980s, shifting these forests from conservation plantings to rotational pine production (Rhodes 2001). NZ's pine forest industry covers >1.7 million hectares (6.3% of total area; 155,600 ha in the Gisborne Region) and contributes ~$3.55 billion to the New Zealand economy (Forest Owners Association & Ministry for Primary Industries 2019).
The forestry industry has positive and negative societal and environmental impacts. For example, there has been strong employment growth (Nixon et al. 2017) and enhanced value through carbon credits (Akpa et al. 2016;Evison 2008;Jindal et al. 2008) but sediment and debris deposition during storm events damages downstream environments (Grant & Wolff 1991;Landcare Research & Scion 2017;Arnold 2018). There are currently additional environmental and safety regulatory pressures on production forestry (e.g., National Environmental Standard for Production Forestry, National Policy Statement for Freshwater Management, National Policy Statement for Indigenous Biodiversity, Health and Safety at Work Act 2015) that may prove untenable for future wood production in some existing areas of pine forest (Clark 2017;Richards 2017). Some of these areas may be more suitable for retirement and natural regeneration post-pine-harvest. A second pathway to permanent carbon sink and natural regeneration is via unharvested pine stands (Climate Change Commission 2021). This scenario has not been addressed in our work as there is insufficient data on the extent of this practice (Forbes & Norton 2021) in NZ, and the perceived negative impacts of pine seedlings and tree fall (Gibson 2021) require further investigation.
Mānuka-kānuka shrubland, composed of mānuka (Leptospermum scoparium J.R.Forest & G.Forst) and kānuka (Kunzea ericoides (A.Rich) Joy Thomps.) is often the first stage of natural regeneration in the Gisborne Region (Williams 1983;Newsome 1987;Wilson 1994;Funk et al. 2009;Overdyck & Clarkson 2012) providing permanent erosion mitigation and possible revenue from mānuka honey. Mānuka-kānuka shrubland accounts for ~70% of regenerating forests in New Zealand (Ministry for Primary Industries 2017b) and native forests have greater potential for erosion mitigation than production pines as clear-cut forest management creates a window of vulnerability for sediment and debris transportation during storm events after harvest (Phillips et al. 2012;Lambie et al. 2018). This is supported by Aburto et al. (2021) who showed native Nothofagus forest had greater erosion mitigation than Pinus radiata production forests in Chile. The scale of work undertaken by Aburto et al. (2021) is comparable with Bergin et al. (1993;1995) and Marden and Rowan (1993) in New Zealand. While native forest will likely have greater erosion mitigation and cleaner water benefits compared to pine production, retirement of pine estates will decrease income from wood production (Hall et al. 2017) and other avenues of income may need to be sourced such as honey, biodiversity, or carbon.
To facilitate informed decision-making by landowners and regulatory authorities, we assessed the viability of converting erosion-prone land currently planted in pines to natural regeneration using high-resolution erosion susceptibility modelling, biophysical modelling of mānuka-kānuka shrubland, and an ecosystem services cost-benefit model. The costbenefit model included scenarios with varying levels of erosion susceptibility and with and without honey production to assess economic viability of a range of scenarios.

Landslide susceptibility and hillslope connectivity
Erosion susceptibility assessments for shallow landslides, gully erosion, and earthflows were undertaken at a 1:25,000. Erosion susceptibility under existing pine plantation forest was assessed. Although different areas of forest were of different ages, we assumed a uniform vegetation cover as the inventory of landslide scars was mapped on pastoral land. Therefore, we were unable to investigate variation in land cover type, and the assessment was limited to lithology and topographic factors. Selection of conditioning factors was based on an understanding of the geomorphic process being assessed. Lithology, slope, and aspect were selected as erosion-conditioning factors, and all have direct physical process relevance for slope stability. The conditioning factors were assessed using bivariate statistics to determine weights and combined in a landslide-index (Juang et al. 1992;Ruff & Czurda 2008;Sciarra et al. 2017) to characterise landslide susceptibility, with rainfall being the cause of landslides. Conditions under which past landslides occurred were used to predict occurrence of future landslides where similar site conditions prevail (Varnes 1984;Soeters & Van Westen 1996;Aleotti & Chowdhury 1999). Spatial distribution of landslides on Tertiary-aged terrain was collated from Betts et al. (2017) who mapped 3,164 landslide scars on similar terrain in the Manawatu region, and Veld and de Graaf (1990) who compiled an inventory of 576 landslide scars on Arai Matawai and Emerald Hills stations ~20 km southwest of Gisborne City following Cyclone Bola (1988). For the Cretaceous terrain, a landslide inventory from the Eastern Ruahine Ranges was used where we assumed that the greywacke bedrock is representative of the Cretaceous terrain found within the Gisborne District. The influence of conditioning factors in susceptibility indices was weighted as per Betts et al. (2017), where slope classes were defined and weights determined by calculating 'densities', or the prior probability of failure within each slope class, based on the location of mapped landslides: where P Prior = P{S} is the conditional probability of having a landslide S and Npix j (Slide) are the number of pixels with landslides in slope class, and Npix j (Total) are the total number of pixels in slope class j. The prior probabilities are normalised to create weights W j within the range of 0-1: Scar-slope relationships were established using a national 15-m digital elevation model (DEM) derived from contour data applied separately for Tertiary and Cretaceous terrains. The prior probability of landslide for the Tertiary terrain was calculated as the mean of the scar-slope relationships, and for the Cretaceous terrain, a single scar-slope relationship was used. Slope aspect was calculated using the 15-m DEM to create nine aspect classes, and the relationship with aspect was assessed on a pixel basis. The conditional probability for aspect was calculated in the same way as it was for slope. Finally, the two conditioning factors (slope gradient for each lithology type and aspect) were multiplied to calculate the landslide Lambie et al. Page 4 N Z J For Sci. Author manuscript; available in PMC 2022 September 08.

EPA Author Manuscript
EPA Author Manuscript EPA Author Manuscript susceptibility index, which was then classified into five categories -none, low, moderate, high, and very high using equal weights.

Hillslope connectivity
Connectivity between hillslopes and waterways assessed if a landslide would deliver sediment to a waterway (Dymond et al. 2006). The 15-m DEM was used to estimate the likely flow-path of material generated by a landslide, its flow direction, and potential intervening accumulation zones, to determine whether sediment and/or slash could potentially enter a stream network. If the flow path encountered any significant flat land (consecutive pixels below four degrees of slope), the source pixel was tagged as 'nonconnected', as we assumed material to be deposited on the flat terrain before reaching the stream. Otherwise, the pixel was tagged as 'connected'. Connectivity was determined for parcels of land with moderate-very high landslide susceptibility as these land parcels were most likely to be retired from pine forest production and reverted to natural forest that affords longer-term erosion control.

Gully and earthflow susceptibility
Gully inventories mapped from aerial photography flown in 1957 (1:15,000) and 1997 (1:26,000) were combined  to identify actively eroding gullies, and the units (with an average size of 110 ha) within the New Zealand Land Resource Inventory (NZLRI) (Landcare Research 2010) that would be most likely to develop gullies. The combined mapped gully layers were then intersected with the NZLRI units. A matrix was used to create gully susceptibility classes (1-5) based on the number of existing gullies (frequency), the potential for future gullies to develop within each NZLRI unit, and the proportion (% of area) of each unit affected by current and past gullying (magnitude). Of the 2,097 NZLRI units in the Gisborne Region identified in the 1970s as having 'present gully erosion', 126 units showed no evidence of gullies when mapped in 1957 and 1997 ). On the assumption that gullies were present in these units in the 1970s, the gully severity ranking assigned to these units at the time of the original mapping of LRI units was adopted.
Earthflow susceptibility assessment draws on the NZLRI (Landcare Research 2010), which includes an erosion severity ranking for earthflows as mapped in 1990 and assigned to classifications between 'none' to 'very high'.

Biophysical modelling
Provision of erosion mitigation by forests is strongly linked to its growth characteristics, which are influenced by climatic, soil and tree species factors. Natural regeneration in New Zealand often begins with mixed mānuka (Leptospermum scoparium and kānuka (Kunzea ericoides var. ericoides) stands (Stephens et al. 2005). Although natural regeneration can begin with broadleaf/podocarp species (e.g., Cameron 1960), there is no published information on the growth and structure of early growth podocarp forest. We, therefore, used mānuka-kānuka shrubland as a representative of natural regeneration as this is the most likely pathway in the Gisborne Region (Funk et al. 2009; Ministry for Primary Industries 2017b).
Growth simulations of mānuka for the Gisborne Region were assessed using a comprehensive dataset of tree growth parameters (e.g., height and diameter at breast height) from across New Zealand from Payton et al. (2010), the New Zealand National Vegetation Survey Databank (NVS) and other unpublished datasets.
Carbon (C) accumulation in mānuka-kānuka shrubland was estimated using model simulations with the physiological model CenW Version 5.0 (Kirschbaum 1999a). CenW has been used extensively to predict the growth of pine stands (e.g., Kirschbaum 1999a, b;Kirschbaum & Watt 2011;Kirschbaum et al. 2012). The model and its source code are available at: http://www.kirschbaum.id.au/Welcome_Page.htm, with a list of relevant equations available at http://www.kirschbaum.id.au/CenW_equations.pdf. For the present work, it was parameterised against the observations from Payton et al. (2010) to simulate the growth of mānuka based on external environmental drivers (temperature and rainfall), stand-internal factors (stand density), and weed competition.
In CenW, weeds are modelled as a separate entity that competes with the main species of interest for nutrients, water, and radiation. Normally, the tree canopy eventually overtops weeds and out shades them, with the time course depending on the initial biomass of weeds and overstorey species and the parameters that define the competitive properties of the weed layer. Here, we assumed a maximal height of the weed layer of 0.5 m and initial weed biomass as 1500 kgDM ha −1 and mānuka-kānuka stands as 200 kgDM ha −1 . These parameters correspond to primary competition by a well-established grass understory.
CenW can be run over periods of many decades with an underlying daily simulation time step. It simulates stand properties and dynamics, such as leaf-area development, stand height, basal area development, litter-fall and exchange of both water and carbon dioxide based on daily inputs of minimum and maximum temperature, solar radiation, rainfall, and vapour pressure (Kirschbaum 1999a). It also requires estimates of site fertility, soil water-holding capacity, and silt and sand fractions as a measure of soil texture. We used 20 years of daily weather input data from the Virtual Climate Station Network (VCSN; National Institute of Water and Atmospheric Research Ltd, Appendix 1). Daily VCSN data were estimated on a 0.05° latitude/longitude grid (Tait et al. 2006;Tait 2008;Tait & Liley 2009) as described by Kirschbaum and Watt (2011). Soil water-holding capacity and the percentage of silt plus clay were obtained from the National Soils Database (Landcare Research 2020). Characteristic climate variables for the Gisborne Region are presented in Appendix Figure A1.
Predictions of different measures of growth were fitted to independent observations from 69 stands located at 52 distinct locations throughout New Zealand (Payton et al. 2010), with growth measures mostly consisting of total stand biomass at respective ages inferred from ring counts. We also used the observed distribution of mānuka-kānuka throughout New Zealand to derive functions to constrain the environmental performance of mānuka-kānuka under more extreme climatic conditions than the sites sampled by Payton et al. (2010).
Some parameter estimates were based on the earlier modelling work of Whitehead et al. (2004) and Whitehead and Walcroft (2005)  Burrows 2 . Other specific information on stand allometric properties was sourced from Scott et al. (2000) and leaf nitrogen concentrations from Ross et al. (2009). The range of parameter values was further constrained to remain within physiologically plausible bounds to retain the physiological integrity of the simulations. The modelling procedure was similar to that described by Kirschbaum and Watt (2011) for pine growth.

Cost-benefit analysis
A cost-benefit analysis using a combination of a spatially explicit agri-environmental economic land-use model (NZFARM) (Daigneault et al. 2018) and other non-market valuation methods was used to monetise changes in land use based on categories of erosion susceptibility. Landslide susceptibility categories 4 (moderate susceptibility with high waterway connectivity), 6 (high susceptibility with high waterway connectivity), 8 (very high susceptibility with high waterway connectivity), and land highly prone to gullying (as per Spiekermann and Marden 2018) were the focus of model scenarios. The model scenarios had two levels of landsliding susceptibility with more categories included in Scenario 2 and consistent gully susceptibility across all scenarios. The erosion susceptibility scenarios were tested with and without honey production as a potential offset to losses in wood profits. Scenario 1, landslide classes 6+8 (high and very high landslide susceptibility with high waterway connectivity) and high gully category converted from pines to native forest with honey production. Scenario 2, landslide classes 4+6+8 (moderate, high, and very high landslide susceptibility and high waterway connectivity) and high gully category converted from pines to native forest with honey production. Scenario 3, same as Scenario 1, without honey production and Scenario 4, same as Scenario 2, without honey production.
Model scenarios were compared to the 'Baseline' scenario, where all identified land remains in pine plantation forestry. The baseline was established using a land use map of the Gisborne Region (AgriBase; Asure Quality 2020) and the New Zealand Land Cover Database (Landcare Research 2015). The model includes assessments of production impacts including profit from land converted from pine forestry, profit from land converted to mānuka honey or related production, value of the land, and planting/native afforestation costs and environmental impacts including carbon cycling and water quality.
The analysis was conducted over 62 years to reflect two cycles of 30-year pine rotation (Hockey & Page unpublished 3 ; Ministry for Forestry 1994). Since the benefits and costs occur over different time periods, the net present value was determined at two discount rates (4 and 6%; Te Tai Ōhanga 2018) to calculate overall impacts: where t is each year (up to 62), and r is the discount rate.
Treasury New Zealand recommended the use 4 and 6% rates to convert future values to present values as money in the future is worth less than the same amount in the present (Te Tai Ōhanga 2018). The range of discount rates also reflects the different objectives of landowners. For instance, the lower end of the range might reflect the objective of achieving intergenerational equity as some landowners (e.g., Māori) might not be driven solely by profits but by other societal goals such as intergenerational resource sustainability as well as maintaining cultural and spiritual benefits for future generations (Goulder & Williams 2012; Harmsworth and Awatere 2013).
Net present value of economic returns for the pine forestry sector were calculated as follows: where p is the log price (Te Uru Rākau 2020), Vol is the log volume harvested in each year (from CenW; Kirschbaum & Watt 2011), and Cost is the fixed and variable costs of production. Fixed and variable costs including logging, cartage, roading, and growing costs were derived from the literature (Olssen et al. 2012). All future values of the production function were then discounted to get total present values. We finally derived the equal annual equivalent (EAE) from the estimated net present value to compare the annual economic returns from pine forestry with other land uses. Economic returns from mānuka honey produced from naturally regenerating mānuka-kānuka shrublands were based on Daigneault et al. (2015).
Environmental factors (E i ) including nutrient leaching, erosion, and carbon sequestration were monetised following estimation on a per hectare basis (γ env ), as impacted by soil type, land cover, and land use. By aggregating the per hectare values of these parameters across the land area under pine (X), we could estimate the total environmental factor outputs from pine forestry: Surficial and stream bank erosion were simulated using the NZEEM model (Dymond et al. 2010), and nutrient values obtained (e.g., Parfitt et al. 1997;Lilburne et al. 2010). GHG emissions were sourced from MPI carbon look-up tables (Ministry for Primary Industries 2017a). To reflect the impact of land conversion from pine production to native forest on the environmental outputs, Equation 5 was modified to: where Z is the area of the land converted to native forest. The parameter γ env specifies the environmental impacts of pine forest after accounting for land conversion, while μ env describes the impact of native forest on the environmental factors. Several non-market techniques have been developed to place a price on changes in water quality, including hedonic pricing (Boyle et al. 1999), recreation demand (Massey et al. 2006), and stated preferences (Moore et al. 2018). Of those techniques, estimates from stated preference studies capture the widest range of people and values, both use and non-use, in their application depending on the scope of the study. We therefore focus on stated preference values and used willingness to pay (WTP) for improved water quality. To monetise water quality improvements associated with each scenario, a benefit transfer for people's WTP for reductions in nutrient leaching was used, which is common technique for determining the value of an ecosystem service (e.g., Aguilar et al. 2018;Huber and Finger 2020;Tian et al. 2020  ) and $6.50 t −1 , where the upper value includes only avoided flood damage and water treatment costs (Barry et al. 2014).
Annual honey profits in mānuka-kānuka shrublands vary largely depending on honey quality ranging between $98 ha −1 for low unique mānuka factor (UMF) to $1000 ha −1 for high UMF (Walsh et al. 2017). These profit estimates were based on a production per hive of 30 kg and 35 kg and a price per kg of $26.30 and $40 for the $98 ha −1 and $1000 ha −1 , respectively (Burke 2015;Wetere 2015). The variation in profits was due to differences in mānuka productivity, and capital and operating costs. UMF is a rating system designed to ensure purity, quality, and authenticity of mānuka honey, and the higher the UMF rating the higher the monetary value of the honey. While no spatially explicit data is available on where high UMF honey has been produced, recent research predicted that soil quality, rainfall, climate, and genotype are critical variables in determining the UMF level in a particular area. To identify areas suitable for high UMF honey, we therefore used temperature and precipitation prediction equations  to define the probability of occurrence of mānuka-kānuka shrubland as per Walsh et al. (2017). The probabilities were summed and assigned to NZFARM polygons based on an area-weighted average to create an index describing the relative likelihood of occurrence of mānukakānuka. Based on conversations with local farmers and associations, the top 10% of polygons was assumed to be areas suitable for high UMF honey and therefore allocated a profit of $1000 per hectare, while the rest of the polygons were assumed to produce low UMF and were allocated a profit of $98 per hectare.
Much of the value of agricultural land is tied to the profits associated with it, so we assume that the opportunity cost of removing land from production is reflected in the change in the NPV of future profits. NZFARM outputs change in net revenue represent Lambie et al. Page 9 N Z J For Sci. Author manuscript; available in PMC 2022 September 08.

EPA Author Manuscript
EPA Author Manuscript the central opportunity costs of the change from pine production to native forest with and without honey production. We were able to monetise several notable impacts associated with changes in land use. However, it is important to note that there are a number of other impacts that we were not able to quantify or monetise (Table 1). Biodiversity was quantified as per Walsh et al. (2019) and presented as 'restored significance' as per Mason et al. (2012) and Carswell et al. (2015), but were not monetised and therefore not included in the economic analysis. Further, the impacts of weed control and competition were not able to be included in the analysis due to limited support information.

Erosion susceptibility
Incidences of landslides increased markedly on slopes >16° in Tertiary terrain and often delivered sediment and forest slash into the nearest watercourse. In contrast, on the Cretaceous terrain, shallow landslides are predominantly restricted to the steep flanks of the Raukumara Range. Slope aspect was an important factor for landslide susceptibility in our study area, as found by others (e.g., Yalcin & Bulut 2007;Galli et al. 2008;Ruff & Czurda 2008;van Westen et al. 2008). We found a disproportionate number of landslides on north to north-east facing slopes, and landslide density on north facing slopes was 50% higher than average.
In the Gisborne Region, 73.2% (545,700 ha) had very low-low landslide susceptibility and 21.1% (157,600 ha) had moderate-very high susceptibility. Moderate-very high areas had landslides that were likely to deliver sediment and other material into nearby waterways. Of the two terrain types, Tertiary underlies the largest proportion of the Gisborne Region and landslide susceptibility was very low-low on 67.4% of hill country areas, and moderatevery high on the remaining 32.6% of hill country, of which approximately 165,000 ha (22%) of hill country slopes are directly connected to waterways (Fig. 2). Overall, 125,000 ha (24.7%) of hill country slopes are moderate-very highly susceptible and have potential for landslides and/or anthropogenic disturbances resulting in sediment, and any associated woody debris, entering a water course.
In comparison, Cretaceous terrain has slopes that were generally less steep and landslide susceptibility was very low-low on 85.5% of hill country areas, and moderate-very high on 14.5% (34,800 ha) of remaining hill country areas, of which 13.7% (32,700 ha) of hill country slopes were directly connected to waterways (Fig. 2).
Land under pine forests has higher levels of landslide susceptibility than the average for the whole region (Fig. 2) as they were established on erosion prone land to mitigate erosion. A larger proportion of the current forest estate is on Cretaceous than Tertiary terrain (Landcare Research 2015). However, land occupied by pine forests on Tertiary terrain has a much higher proportion of land susceptible to landslides (35.9% in classes moderate-very high) than the pine forests on Cretaceous terrain (26.8%).
Within the Gisborne Region, 19% of the region is moderate-very high for gully erosion susceptibility (Fig. 3a). 22.5% of hill country is classified as susceptible to earthflow

EPA Author Manuscript
EPA Author Manuscript erosion, with just 8.8% classed as moderate-very high susceptibility (Fig. 3b) and is associated with areas dominated by mudstone or crushed argillite.

Tree growth
To better understand the growth potential of mānuka-kānuka stands, we parameterised the CenW model against the observed data of Payton et al. (2010) (Fig. 4). Stands showed some moderate growth potential over the first 20 years, reaching stand biomass of about 50 tC ha −1 . Growth rates then slow, with peak stand biomass typically reached by about 50 years. Stands then tend to degenerate over further time or are invaded by taller trees that then dominate. The observed growth patterns could be captured well by the CenW simulations (Fig. 4b), with a calculated model efficiency of 0.65.
Growth potential of stands is strongly affected by initial stand density, which interacts with weed competition. Once seedlings over-top the competing grass layer, a period of rapid growth can commence. Stand growth rate tends to be reduced through self-thinning that reduces the number of living trees with associated loss of stand biomass. Growth rates are therefore typically only a fraction of those achieved by well-managed P. radiata stands (cf. Kirschbaum et al. 2011).
Analysis of growth rates (Payton et al. 2010) and the patterns of temperature and rainfall in relation to the natural distribution of mānuka-kānuka stands across the country were used to estimate tree growth in response to environmental drivers. We had insufficient information to make a distinction between the growth of mānuka and kānuka and assumed the same environmental limitations for mānuka stands and mānuka-kānuka shrublands. Mānuka-kānuka stands require mean annual temperatures of more than ~5°C (Appendix Fig. A2a), and growth increases rapidly with increasing temperature to about ~12°C (Fig.  A2a). Similarly, mānuka-kānuka stands require annual rainfall > ~300 mm yr −1 (Appendix Fig. A2b), and growth increases steeply for optimal performance at ~800-1000 mm yr −1 , with no further growth response at higher rainfall. However, it is not currently known whether extremely high rainfall beyond 2,000 mm yr −1 adversely affects the growth and persistence of stands.
For the Gisborne Region, simulated growth of mānuka-kānuka was highest in the north-east, where it could reach up to 1.2 tC ha −1 yr −1 , and marginally lower in the coastal regions on both the north and southeast (Fig. 5a). Growth was slightly lower along the higher-elevation inland ridge where cooler temperatures limited growth to 0.7-0.8 tC ha −1 yr −1 . Low values were also obtained for a few southern coastal sites where low precipitation combined with soils with low soil water-holding capacity to cause water-stress limitations.
Pines are capable of much faster biomass and wood accumulation than naturally regenerating shrubland. Stem diameters in naturally regenerating stands are much smaller than for pines because of the much higher stand density of naturally regenerating stands, so that their biomass must be spread over a larger number of stems (Fig. 5). Pine growth was modelled to be highest in the regions with moderately high rainfall, with biomass growth of up to 12 tC ha −1 yr −1 (Fig. 5b). This corresponds to wood production of 350-450 tDM ha −1 over a 30-year rotation (data not shown). In contrast, growth in the high-rainfall Raukumara Range ridge running from north-east to south-west was slightly lower with biomass growth of only 7-8 tC ha −1 yr −1 . Unlike mānuka (Stephens et al. 2005), pines generally do not grow well in regions with excessive rainfall (Kirschbaum & Watt 2011).

Cost-benefit analysis
The estimated area of pine forests for conversion in each scenario was 315 and 556 ha, or 0.33 and 0.59% of the total planted pine area in the study area. In the first two scenarios, which included honey production in native forest areas, 1-2 ha was classified as high-UMF production ( Table 2) and were estimated to have high economic returns, while the remainder of the land was assumed to produce lower-value honey with poorer economic returns. Seven and 15 ha were estimated to incur additional expenses for successful native regeneration (Table 2), including purchasing plants and planting costs. The remaining land within the scenarios was predicted to have sufficient landscape and geographic factors to permit natural regeneration (Walsh et al. 2019) and was likely to be mānuka-kānuka shrubland initially (Bray et al. 1999;Funk et al. 2009; Ministry for Primary Industries 2017b). There was lower reduced carbon storage and erosion under all four scenarios compared to the baseline (Table  2). There were also some changes in biodiversity and water quality but given the size of those changes and the lack of usable nonmarket values, they were difficult to monetise. Restored significance, as an estimate of biodiversity improvements, was 324±20 ppb for all scenarios.
Overall, the monetised NPV of converting from pine in each of our scenarios is negative (Table 3). Converting pine forestry land reduces profits from production and the land value, and the two combined can represent significant decreases. Although revenue can be obtained from other enterprises based on natural afforestation, such as mānuka honey or oil, these are not as profitable as growing pines. However, there were some environmental benefits that could not be monetised and therefore were omitted from the NPV calculations. Highly erosion-prone land is usually on steep slopes and the analysis of profits and land values may be overestimated (Table 3).
Across the four scenarios, losses in NPV ranged between $3 and $8 million over the 62-year period or ~$150-$250 ha −1 yr −1 and lost profits from forestry represented the largest costs of conversion. Decreases in carbon sequestration were also notable, as native forest stores less carbon across 62 years than pines. Other social benefits, such as environmental and cultural factors, are more difficult to quantify. Yao et al. (2014) found that there is a willingness of the New Zealand public to pay for enhancement of forest ecosystems for biodiversity provision, but further work is needed to adequately provide values for biodiversity and cultural parameters, which may further offset wood derived profits under our scenarios. Several other categories could be neither quantified nor monetised. In terms or overall impacts, these omitted categories would likely decrease the overall negative impacts of our scenarios. These omitted values may provide compelling reasons to promote native afforestation. Cultural impacts, for instance may justify the loss in forestry profit is some cases.
It is also difficult to calculate potential changes in employment resulting from our modelled land use conversions. However, since the scenarios we analysed affect only small land areas, the employment impact should be minor. Larger conversions would likely result in proportionately larger employment changes with wider regional economic impacts.

Discussion
Both plantation pine and natural regeneration have long-term benefits for catchments and communities. However, clear-felling of pine stands can lead to introduction of sediment and woody debris to waterways during storm events (Marden et al. 2006;Imaizumi et al. 2008;Payn et al. 2015;Issaka & Ashraf 2017;Spiekermann & Marden 2018). Erosion susceptibility is inherently variable within catchments and geological terrains suggesting that land cover should ideally be as diverse as the landforms to create a mosaic of land uses across a property or catchment.

Economic impacts of conversion from pine production to native forest
If assessing the viability of converting pines to natural regeneration on a purely economic basis, pines are considerably more profitable than natural regeneration, even when other market benefits of natural regeneration are included. The costs of conversion would be borne by pine forest owners, in terms of lost profit and changes in land value. Conversely, many of the benefits from forests accrue to the general population, such as improvements in local air and downstream water quality, carbon sequestration, and cultural impacts (e.g., Nowak et al. 2012).
Mānuka-kānuka shrubland is often the first phase of natural regeneration (Overdyck & Clarkson 2012), and under suitable circumstances, honey profits may be used to offset some of the losses from retiring land under pines. Scenario 2 may represent the ideal situation whereby the most highly vulnerable parts of the landscape (susceptible to both landslides and gullying) are retired from production forestry and allowed instead to develop a permanent tree cover, with the associated benefit of honey production.
Many historical pine plantations were established in very difficult, steep terrain, often far from urban centres, ports, and mills increasing operational and infrastructural costs (Raymond 2012). These impacts on profits were not included in our analysis and strategies are being developed to improve harvesting efficiency on steep sites, it is likely that profits for these slopes may be underestimated (Raymond 2012(Raymond , 2014Amishev et al. 2013). On particularly difficult slopes, trees may remain unharvested, especially during periods of low log prices (Review Panel 2014) and therefore act as a nursery crop of natural regeneration and any profits for these sites will be unrealised.
Weeds can reduce the growth and survival of target species and is an important factor for long-term erosion mitigation and monetary profits. Weeds can initially compete strongly with mānuka-kānuka seedlings for light, water, and nutrients, but once mānuka trees are taller than competing weeds, they gain largely unrestricted access to solar radiation and will out-shade low-stature weeds once the canopy is sufficiently large and dense. Our growth simulations were based on assuming a well-established grass layer as the principal weed competitor. Because of its high initial biomass, grass acts as an effective competitor in the early establishment phase, but because of its low maximum height, mānuka-kānuka

EPA Author Manuscript
EPA Author Manuscript seedlings could soon out-shade the grass layer and become the dominant plant type.
If shrubs, such as gorse or broom, were the principal competitive weed species, their competitive inhibition of the growth of mānuka-kānuka seedlings could have persisted for longer or could even have prevented the establishment of the shrubland altogether.
These impacts of weeds have been well characterised for pines (e.g., Richardson et al. 1996;Watt et al. 2007), but we are not aware of systematic studies of weed effects on the establishment of mānuka stands. Pine seedlings are another substantial weed problem occurring post-harvest and during the natural regeneration stage. Pine seedlings can have a substantial impact on native species establishment and survival (Marlborough District Council et al. 2016). Further, the costs of pine seedling control can be substantial ($150-500 ha −1 ) and affect the species composition of naturally regenerating forests as different native species differ in their sensitivity to wilding control with herbicide applications (Marlborough District Council et al. 2016;.

Drivers for conversion of pine production to native forest
Land-management decisions are also subject to social influences as well as the economic factors (Miller et al. 2007;Bhandari et al. 2015). Enhancing biodiversity values is a driver for shifts from pine to native forests (Marlborough District Council et al. 2016;. Pine forests have good biodiversity values, providing habitat for many plants and animals including threatened species (Brockerhoff et al. 2001;Bremer et al. 2010;Pawson et al. 2010;Michelsen et al. 2014;Berndt & Brockerhoff 2019). Native forests typically contain a greater species richness compared with pine forests, as native forest has a greater range of food sources that are unimpacted by harvest and forest structural diversity providing a greater range of habitat; however, native forest is highly temporally and spatially variable (Clout & Gaze 1984;Díaz et al. 2005;Pawson et al. 2010). As a result, some landowners who highly value these attributes may be willing to reduce their profits from pine revenues.
Biodiversity offsetting is also a potential pathway for recognising the value associated with native forests (Department of Conservation 2014). However, a monetary or currency value associated with biodiversity offsetting in New Zealand remains highly complex and a valuation framework is required (Department of Conservation 2014). Biodiversity offsetting is generally supported through resource management undertaken by Councils, under the umbrella of the Resource Management Act. Substantial work has been undertaken to establish guidance on biodiversity offsetting (Department of Conservation 2014; Maseyk et al. 2018) and it is possible that this mechanism will be further supported under the Resource Management Act reforms and in the future (Ministry for the Environment 2021).
There are also pathways where financial value can be assigned to projects that enhance biodiversity. For example, Green bonds, Green funds, sustainability linked loans and biodiversity credits (Chartres 2021

EPA Author Manuscript
EPA Author Manuscript require significant development of a framework in which to undertake credit trading and recognition in New Zealand.
Shifts from pine to native forest cover can also be driven by cultural values (Hēnare 2014). Former state-owned forest assets are returning to Māori under Treaty of Waitangi claims and could result in 41% of post-harvest pine forests being Māori-owned (Miller et al. 2007). Māori connect with indigenous forests on a spiritual level associated with whakapapa (genealogy) and kaitiakitanga (guardianship), which provide additional noneconomic factors for decisions on land use (Miller et al. 2007) that may be attractive to Māori landowners. In the Gisborne Region much of the most erosion-prone land is on Māori-owned land (Miller et al. 2007), which may provide an added incentive for increasing indigenous forests.
Carbon credits may also drive natural regeneration. However, changing from a fast-growing tree crop such as pine trees to potentially slower growing native trees will reduce carbon stocks and carbon sequestration and ultimately result in less carbon credits compared to pine trees (Kimberley et al. 2014) and is more profitable when converting from pasture to forests. Further, carbon accrual in regenerating forests is calculated generically from look-up tables which are particularly lacking with respect to native species, and do not contain regional or species-specific information (Ministry for Primary Industries 2017b). Carver and Kerr (2017) suggest that updating the tables to include more native tree specific information will facilitate inclusion of native forests in the ETS and therefore greater recognition of these forests for carbon income.

Conclusions
We suggest a mosaic of land use within a property (or catchment) may be the best overall option, where moderate to highly susceptible land is left to naturally regenerate post-pineharvest, while pine production is maintained on land less susceptible to erosion. For this multi-land use approach to be considered viable and adopted successfully, landowners will need to understand the erosion susceptibility of their land to create an impetus to shift from pine trees to native forest. This would benefit from support from local and central government agencies with the overall aim to meet cleaner water aspirations. Many of the current mechanisms to encourage erosion mitigation are not applicable to those wanting to shift from pines to native forests and only support planting in previously unforested areas (e.g., One Billion Trees Programme, Erosion Control Funding Programme). Further, the relatively lower amount of carbon credits that native forests accrue relative to pine plantations may also discourage shifting to native species. Overall, these results illustrate the economic value of the pine forestry industry, and the consequent challenges with transitioning to other land uses despite the greater environmental benefits associated with reducing erosion.

Gisborne climate
Solar radiation across the Gisborne region is fairly uniform, with just slightly higher values at the coast and slightly lower values inland at higher elevations (Fig. A1a). Temperature shows greater variation, with mean annual temperatures above 15°C at the north-eastern coastal tip and dropping to below 10°C in the higher inland locations, especially towards the south-west (Fig. A1b, c, d). There is rainfall of 2,500 to over 3,000 mm yr −1 in the upland regions, but less than 1,500 mm yr −1 in most of the coastal regions, and rainfall of only around 1,000 mm yr −1 inland from Gisborne city (Fig. S1e). In terms of temperature extremes, the area near Gisborne city has experienced the greatest temperature extremes of recorded maxima above 35°C, while the northern side of the peninsula has reached only 26-28°C. Absolute minima have generally been lowest in the highland region with values below −7°C, whereas the coastal regions generally remained milder with minima mostly not falling below −3°C, except for a large part of the east coast where temperatures have fallen to below −5°C.

FIGURE A1:
Key environmental variables and soil characteristics of the Gisborne Region, radiation (a); mean temperature (b); recorded maximum (c); and minimum temperatures over a 20-year period (d); rainfall (e); and soil water-holding capacity (f). Climatic data refer to the 1980-1999 period, with panels giving either the averages recorded over that period (a, b, e) or the absolute extremes recorded over that period (c, d). All data are shown at 0.05-degree resolution.

FIGURE A2:
The relative growth response of mānuka-kānuka stands for different mean annual temperatures (a); and mean annual rainfall (b). These relationships have been developed based on analysis of the data of Payton et al. (2010) and the observed distribution of mānuka-kānuka stands across New Zealand.

List of abbreviations
UMF unique mānuka factor NPV net present value   Landslide susceptibility, as a percentage of: (a) major catchments and Gisborne Region; and (b) for the area of Tertiary and Cretaceous terrains, including pine forest areas only. For landslide susceptibility, the percentage of hill country with moderate to very high susceptibility and potential to deliver sediment (connectivity) to waterways is also shown.    Observed and modelled growth potential of mānuka-kānuka stands of different ages, showing observed (symbols) and modelled (solid curve) biomass (a) and observed biomass plotted against modelled data (b). Observed data are from Payton et al. (2010) and simulations from CenW.  Simulated growth over 20 years of: (a) total biomass for naturally regenerating mānukakānuka stands (100,000 stems ha -1 ); and over 30 years of: (b) total biomass for pine (P. radiata) in the Gisborne Region. Note the 10-fold difference in scales between mānukakānuka stands and pine.   Model inputs for cost-benefit valuation for four scenarios, two scenarios with high (6) to very high (8) landslide susceptibility and two with moderate (4), high (6), and very high (8) landslide susceptibility, all scenarios included high gully prone land. Scenarios 1 and 3 are the same except for the inclusion of honey in Scenario 1. Scenarios 2 and 4 are the same except for the inclusion of honey in Scenario 2. LS is the landslide susceptibility category, and Δ indicates the change from the baseline scenario.  Baseline NPV (NZ $) across 62 years and change in NPV for four scenarios (4% and 6% discount rates), two scenarios with high (6) to very high (8) landslide susceptibility and two with moderate (4), high (6), and very high (8) landslide susceptibility, all scenarios included high gully prone land. Scenarios 1 and 3 are the same except for the inclusion of honey in Scenario 1. Scenarios 2 and 4 are the same except for the inclusion of honey in Scenario 2. Where LS is the landslide susceptibility category, 'seq' represents sequestration, 'afforest' represents afforestation, and Δ indicates the change from the baseline scenario.