Natural forests in New Zealand–a large terrestrial carbon pool in a national state of equilibrium

Thomas Paul,Mark O.Kimberley and Peter N.Beets

Abstract Background:Natural forests cover approximately 29%of New Zealand’s landmass and represent a large terrestrial carbon pool.In 2002 New Zealand implemented its first representative plot-based natural forest inventory to assess carbon stocks and stock changes in these mostly undisturbed old-growth forests.Although previous studies have provided estimates of biomass or carbon stocks,these were either not fully representative or lacked data from important pools such as dead wood (coarse woody debris).The current analysis provides the most complete estimates of carbon stocks and stock changes in natural forests in New Zealand.Results: We present estimates of per hectare carbon stocks and stock changes in live and dead organic matter pools excluding soil carbon based on the first two measurement cycles of the New Zealand Natural Forest Inventory carried out from 2002 to 2014.These show that New Zealand’s natural forests are in balance and are neither a carbon source nor a carbon sink.The average total carbon stock was 227.0±14.4 tC·ha?1(95%C.I.) and did not change significantly in the 7.7 years between measurements with the net annual change estimated to be 0.03±0.18 tC·ha?1·yr?1.There was a wide variation in carbon stocks between forest groups.Regenerating forest had an averaged carbon stock of only 53.6±9.4 tC·ha?1 but had a significant sequestration rate of 0.63±0.25 tC·ha?1·yr?1,while tall forest had an average carbon stock of 252.4±15.5 tC·ha?1,but its sequestration rate did not differ significantly from zero (?0.06±0.20 tC·ha?1·yr?1).The forest alliance with the largest average carbon stock in above and below ground live and dead organic matter pools was silver beech-red beech-kamahi forest carrying 360.5±34.6 tC·ha?1.Dead wood and litter comprised 27% of the total carbon stock.Conclusions:New Zealand’s Natural Forest Inventory provides estimates of carbon stocks including estimates for difficult to measure pools such as dead wood and roots.It also provides estimates of uncertainties including effects of model prediction error and sampling variation between plots.Importantly it shows that on a national level New Zealand’s natural forests are in balance.Nevertheless,this is a nationally important carbon pool that requires continuous monitoring to identify potential negative or positive changes.

Keywords:Forest inventory,Natural forests,Carbon stocks and change,Terrestrial carbon pool

Introduction

Natural temperate forests contain a significant proportion of the non-atmospheric land-based carbon,with mature primary forests in the southern hemisphere forming some of its most carbon dense examples (Keith et al.2009).Maintaining these significant carbon stocks at high levels is seen as paramount for mitigating climate change (Carey et al.2001;McKinley et al.2011;Keith et al.2014).Their reduction either through land-use change,unsustainable forest management or human induced degradation could result in large emissions that would be difficult to recoup (Keith et al.2014;Lewis et al.2019).It is evident that the accurate estimation of carbon stocks in natural forests is important for countries with large areas of natural forests lacking direct human impacts.This carbon pool is potentially at risk from indirect or direct human impacts,and from natural disturbances (Kurz et al.2008;Li et al.2011;Brando et al.2014).Country-level estimates can be made using national forest inventories which can provide critical insights concerning carbon stock changes,and when synthesised together can provide global change estimates(Keith et al.2009).

New Zealand is in the cool temperate moist zone,where high carbon stocks are expected (Lal and Lorenz 2012).Forest covered as much as 85% of its land area prior to human habitation (McGlone 1989) but has been greatly reduced by fire and forest clearance since the arrival of humans approximately 800 years ago (Beaglehole 2012).Currently,natural forest monitored by the forest inventory described in this paper,covers 29% of New Zealand’s land area.It consists mainly of primary or oldgrowth forests but also includes secondary earlysuccessional forests and shrublands growing on land that was historically cleared for farming or burned,and which can be distinguished by species composition,being characterised by seral shrub and tree species (Wiser et al.2011;Holdaway et al.2017).Prior to 1987,it was managed by the New Zealand Forest Service for either timber production,utilizing desirable tree species,or soil conservation purposes and therefore protected.When the New Zealand Forest Service was disestablished in 1987,most harvesting of natural forest in New Zealand had ceased,and the Department of Conservation (New Zealand Government 1977) was established to manage the forest for non-timber values.

To evaluate the role of New Zealand’s natural forest globally in the mitigation of climate change through carbon storage and sequestration,and to meet its international reporting commitments under the United Nations Framework Convention on Climate Change(UNFCCC),the Kyoto Protocol and more recently the Paris Agreement,New Zealand began implementing its first nationally representative natural forest inventory in 2002 as part of the wider Land Use and Carbon Analysis System (LUCAS).The natural forest inventory covers pre-1990 natural forest which is defined as land that was covered in naturally occurring (i.e.,not planted) forest in 1990,and remains in forest.Separate forest inventories cover New Zealand’s post-1989 natural forest (Beets et al.2014),pre-1990 planted forest (Beets et al.2012a),and post-1989 planted forest(Beets et al.2011).The natural forest inventory is designed as a representative network of permanent plots laid out on an 8 km ×8 km grid covering the mapped area of pre-1990 natural forest and allows the construction of unbiased estimators of carbon stocks and biodiversity metrics with known precision (Coomes et al.2002).Most of the network plots were newly installed as no systematic national inventory of natural forest existed.However,18% consisted of previously installed permanent plots located close to grid points.To ensure consistency the newly installed plots used the same plot installation methods as the existing plots.However,a horizontal 20 m radius plot was added to the design to measure carbon in larger trees.The historically installed 20 m×20 m square plots were not slope corrected.Despite the potential limitations of this approach,it was decided to utilize plots with previous measurements within the developing natural forest inventory (started in 2002) to provide estimates for other purposes including the 1990 baseline of carbon stocks(Hall et al.2001).The wall-to-wall satellite-based map used to define the area covered by the inventory is periodically updated and corrected.The analysis described in this paper is based on the map as it existed in 2012(Ministry for the Environment 2020).

Carbon estimates for New Zealand’s temperate natural forest were first provided by Hall et al.(2001) who reported carbon in above ground biomass (AGB) and below ground biomass (BGB) in live trees and shrubs to average 180 tC·ha?1utilising data from numerous National Vegetation Survey plots (Wiser et al.2001).Repeated measurement from these plots showed fluctuations in above ground live biomass at decadal time intervals.However,this survey did not include measurements of dead organic matter pools,and moreover,was focussed around selected catchments (Allen 1993),and therefore lacked the full spatial coverage of natural forest in New Zealand.An analysis of the first completed measurement cycle of the natural forest inventory carried out between 2002 and 2007 was reported by Beets et al.(2009).This analysis found that the total carbon stock in pre-1990 natural forest and shrubland averaged 173 tC·ha?1.The stock in natural forest excluding shrubland was 218 tC·ha?1.These results indicate that New Zealand’s natural forest are significant,holding more carbon per hectare than boreal forests and similar levels to tropical forests (Keith et al.2009).

Analysis of stock change based on the natural forest inventory became possible after the completion of the second measurement cycle in 2014.An analysis of biomass stocks and stock changes from 874 grid locations with two measurements was reported by Holdaway et al.(2017).They found that 16%of inventory plots located in secondary forest showed a significant positive change in total biomass (live above and below-ground biomass,dead wood and litter) with biomass increasing by 2.78 t·ha?1·yr?1.However,the 84% of plots in old-growth forest had an average change of only 0.28 t·ha?1·yr?1.Across all forest,the net change in biomass was estimated to be 0.67 t·ha?1·yr?1and did not differ significantly from zero.This suggested that New Zealand’s natural forest is in a steady state with respect to biomass cycling.

This finding contrasts with results from national forest inventories globally which suggest that forest biomass stocks in forests are generally increasing (Carey et al.2001;Luyssaert et al.2008;Pan et al.2011).In fact,sequestration of carbon by the world’s forests is currently believed to be a major contributor to the global residual terrestrial carbon sink,calculated as the difference between all known carbon sources and sinks (Houghton et al.2018;Ciais et al.2019;Friedlingstein et al.2019).The discovery that intact natural forests are apparently acting as carbon sinks was initially considered surprising as it was generally assumed that mature forests would be carbon neutral (Jarvis 1989;Carey et al.2001) as appears to be the case for New Zealand’s natural forest.Multiple reasons why forests may generally be acting as carbon sinks have been suggested.Models suggest that the effect could be mostly due to CO2fertilisation resulting from rising atmospheric CO2(Schimel et al.2015),with other causes likely to include the effect of N deposition (Du and de Vries 2018),synergies between N and CO2(O'Sullivan et al.2019),the effects of climate change such as warming and lengthened growing season(McGuire et al.2018),and recovering growth in immature forests (Birdsey et al.2006;Pugh et al.2019).Although several of these factors are unlikely to be significant for New Zealand’s forests (e.g.N fertilisation and recovering growth),CO2fertilisation would have been expected to occur in New Zealand as much as in those of other countries and the apparent absence of carbon sequestration in New Zealand’s forest is therefore surprising.We note,however,that the level of sequestration provided by old-growth forests globally is the subject of ongoing debate and recent analysis has cast doubt on its extent (Gundersen et al.2021).While we do not have a good understanding or measured evidence of how these global drivers may affect New Zealand natural forest’s ability to sequester carbon,it is of importance to improve our estimates of carbon stocks and changes to be able to detect such effects.

Given the importance of estimating changes in forest carbon stocks,it is essential that such estimates be made with as much precision and accuracy as possible.The estimate of biomass change in New Zealand’s natural forest provided by Holdaway et al.(2017) was based solely on data from the smaller 20 m×20 m inner plots used in the natural forest inventory and did not include measurements from the larger 20 m radius circular plots used to sample larger diameter trees and dead wood,which are an integral part of the nested plot design of the natural forest inventory (Coomes et al.2002).It has subsequently emerged that large stems were over-sampled in the 20 m×20 m plots,meaning that estimates of carbon stocks based solely on them are too high (Paul et al.2019).It is therefore necessary to provide revised estimates utilising data from the 20 m radius plots.

Here we present revised estimates of per hectare carbon stocks and stock changes in above-ground live biomass (AGB) consisting of living stems,branches and foliage,below-ground live biomass (BGB) consisting of living roots,and dead wood greater than 10 cm diameter(or coarse woody debris,CWD) not contained in the litter and either standing,lying on the ground,or in the soil.We also present estimates of carbon stocks but not stock changes in litter (or fine litter) including the litter,fumic,and humic layers.However,we do not present estimates for soil carbon.All estimates are for New Zealand’s pre-1990 natural forest based on data from two completed measurement cycles of New Zealand’s natural forest inventory.Because the natural forest inventory fully represents New Zealand’s natural forest,it avoids the risk of sampling bias found in some previous studies due to selective sampling or the local focus of single or multiple studies (Fisher et al.2008).Our carbon stock estimates exclude the effects of known biases that resulted from the 20 m×20 m plot placement and layout(Paul et al.2019),correct for a systematic underestimation of the dead wood carbon due to known difficulties in measuring dead wood carbon in natural forests(Kimberley et al.2019),and incorporate recent improvements in assumptions concerning below ground carbon estimation (Easdale et al.2019).

Methods

New Zealand’s natural forest inventory

New Zealand’s Natural Forest Inventory is based on a network of permanent plots installed on a randomly placed systematic 8 km ×8 km grid covering the mapped area of natural forest and shrubland (Fig.1).Satellite imagery was used to derive the mapped natural forest area as at January 1990 which remains forested.The map is periodically updated (Ministry for the Environment 2020) but remained essentially unchanged during the first two measurement cycles of the inventory.The 2012 version of the map was used in the current analysis.The mapped area of natural forest covers 7.74×106ha and is intersected by 1195 grid locations which are used as sampling points for the inventory.

Fig.1 Natural forest inventory grid locations plots showing plots in tall forest(measured twice(blue),and once(red)),and regenerating forests and shrublands(measured twice(yellow),once(gold)).Two locations with no measurable stem data could not be classified into forest types(black squares)

For the first measurement cycle undertaken during 2002–2007,1116 points were randomly selected from the grid locations for field measurement,with resource constraints precluding measurement of the remaining 79 points.The selected points therefore formed a random sample of locations within the mapped area.Of these,access was denied to 47 locations,33 were inaccessible,and the remaining 1036 points were measured.A different procedure was adopted to select locations for measurement during the second measurement cycle undertaken during 2009–2014.Plots that had been sampled in the first cycle were classified into two strata,‘Shrub’ and‘Forest’,based on their species composition.Because there was particular interest in the ‘Shrub’ stratum,all but two of the 79 points in this stratum were remeasured.Of the 957 points in the ‘Forest’ stratum,835 were randomly selected for remeasurement.In addition,10 locations from the 79 that had not been sampled during the first cycle were selected randomly for measurement in the second cycle.In all plots selected for measurement in both cycles,measurements of live stems and dead wood were carried out as described below.However,litter measurements were only obtained in a random subsample of 250 plots in the first cycle and no litter measurements were made during the second cycle.

Plot design and carbon related field measurements

A permanent nested plot design was deployed with a 20 m horizontal radius circular subplot,and a 20 m×20 m square subplot nested in the centre of the circular plot and permanently marked for easy relocation (Payton et al.2004a;Ministry for the Environment 2018).Trees within a 20 m horizontal radius were identified by measuring the horizontal distance with a Vertex between their midpoint diameter at breast height (DBH,1.35 m)and the plot centre peg and were permanently tagged.In the case of marginal trees,the distance to the centre point of each side of the stem was checked.The 20 m×20 m plots were installed by laying 20 m tapes along the ground surface without correcting for slope.This nonhorizontal field method was based on the historic approach used in vegetation surveys undertaken in New Zealand,which included also the permanent tagging of individual trees greater than 2.5 cm DBH (Payton et al.2004a).Horizontal distances measured between permanently marked corner-pegs and angles were used to calculate the true horizontal area using the Bretschneider formula for the area of a convex quadrilateral (Paul et al.2019) resulting in an average“square”plot area of 0.0345 ha across the inventory.Large diameter stems and dead wood with a DBH or minimum diameter respectively of 60 cm were measured within the horizontal 20 m radius plots (0.1257 ha).Live stems with a minimum of 2.5 cm DBH and dead wood with a diameter ≥10 cm were measured within the inner square plot only.

Live tree measurements

All live tree stems in the permanent plot that exceeded the DBH limits for the relevant subplot were tagged permanently.The DBH,species,and status (dead or alive)were recorded at each measurement.A subset of up to 15 more-or-less straight (not excessively bowed) trees from each of the following species groups,hardwoods,softwoods,tree ferns,and palms,if present,was selected across the DBH range and measured for total stem length.If stem lean exceeded 10 degrees from vertical,the lean angle was recorded following standard protocols(Payton et al.2008;Ministry for the Environment 2012)to calculate true stem length.

Dead wood measurements

All dead wood exceeding the diameter limits for the relevant subplot and with a minimum length of 30 cm was tagged depending on its state of decay (Payton et al.2008) and DBH and height recorded for standing dead wood and two end diameters and length recorded for fallen dead wood.The species or genus was recorded(when possible) and a decay class was assigned for each individual dead wood piece (Coomes et al.2002).Stems with status“l(fā)ive”at Time 1 and identified as“dead”at Time 2 (standing or fallen) kept their tag for tracking over time and were remeasured and assigned a decay class.

Quality assurance

Training and annual audits are part of the New Zealand Natural Forest Inventory to ensure field procedures described in natural forest inventory-specific field manuals are followed and applied in a consistent manner over the inventory cycles.Field audits were conducted on 10% of plots in the first inventory cycle (2002–2007) and on 5%of plots during the second inventory cycle (2009–2014).The audits identified a number of areas for improvement,and the data collection manual was revised to provide improved clarity over time.

Prior to undertaking carbon stock calculations,the data were thoroughly checked and corrected,ensuring mortality,ingrowth and individual stem growth was accurately accounted for through time in live tagged stems.A preliminary analysis of the dead wood data from both inventory cycles was also undertaken (Kimberley et al.2019) and supported the adjustment approach to the dead wool pool calculations,described below.

Height estimation

The estimation of carbon stocks and their changes was based on metrics calculated at the individual stem and piece level.This required individual tree heights which when combined with DBH measurements were used to estimate volumes of individual stems which were then used to estimate carbon accounting for wood density,and for dead stems,decay.

As heights were only measured for a sample of stems,unmeasured heights were estimated using height/DBH regression models fitted to height measurements of all live,non-leaning stems measured twice (i.e.,in both 2002–2007 and 2009–2014).Height measurements of leaning stems undertaken during the first measurement cycle were considered to be unreliable and were therefore not used.Furthermore,as described in Kimberley and Beets (2016),tests suggested different criteria were used to select height trees in each measurement cycle.Therefore,to ensure that height changes represented genuine growth effects,only non-leaning stems measured in both measurement periods were used to develop the models.For these non-leaning trees,tree height was assumed to be the measured stem length.The method used the following underlying relationship between height,H (m),and DBH (cm):

with the a and b coefficients adjusted to account for plot and species differences.Full details of this methodology are presented in the supplementary material.

Live stem volumes and resulting carbon stocks

Over-bark volumes (m3) of all live stems measured for DBH other than tree ferns,cabbage trees and palms were estimated using the allometric equations in Beets et al.(2012b) from DBH (cm) and predicted height (m).Dry weight of live stems was estimated by multiplying stem volume by whole stem wood plus bark density.Values by species range from 288 to 770 kg·m?3,and are tabulated in the Supplementary Material of Beets et al.(2012b).For species with no tabulated density,the mean density of the genus,or failing that of the growth type(canopy tree,subcanopy tree,tree-fern/cabbage tree or shrub) was used.Stem carbon was then estimated by multiplying the dry weight by the carbon fraction,assumed to be 0.51 for gymnosperms and 0.48 for broadleaf and other species(Table 4.3 in Volume 4 of Intergovernmental Panel on Climate Change(IPCC)2006).Carbon in branches and foliage was estimated directly from DBH using the equations in Beets et al.(2012b).

For tree ferns,cabbage trees and palms,carbon was estimated directly from DBH (cm) and predicted height(m) using the equations in Beets et al.(2012b).The carbon in below-ground biomass in each live tree or shrub was calculated as a ratio (root/shoot ratio) of the carbon in above ground biomass (stem+branch+foliage).The ratios used were 0.234 for angiosperm trees ≥5 cm DBH,and monocots (palms and cabbage trees),0.194 for tree ferns,and 0.245 for gymnosperm trees and shrubs (Easdale et al.2019).

Dead wood volumes and carbon stocks

For dead standing trees,the original unbroken volume(m3) was predicted using the same allometric equation used for live stems based on DBH and predicted height.The fraction of stem volume remaining after accounting for breakage was then predicted from the ratio of measured height to predicted height by integrating the taper function described in Beets et al.(2012b).As dead wood or Coarse Woody Debris is defined as being at least 10 cm in diameter,the function was only applied to stems≥10 cm DBH,and volume was calculated to the lesser of the measured height and the height corresponding to a stem diameter of 10 cm estimated from the taper function.Volume (m3) of dead wood in stumps was estimated from the stump small end diameter (SED,m) and height (H,m) using the formula for a cylinder:

Volume (m3) of fallen dead wood was estimated from large and small end diameters (LED and SED,m) and piece length (L,m) using the formula for a truncated cone:

Dead wood carbon (kg) was estimated based on volume of the individual fallen dead wood pieces,dead standing trees,spars and stumps using the same decay modifiers and carbon fraction:

In this equation,VDead_Wood(m3) is volume of the standing dead stem,stump,or fallen piece,Density(kg·m?3) is whole stem wood plus bark density,as used for live stems,Decay_Modifier is an adjustment factor applied to allow for the loss in whole stem density due to decay and fragmentation.Generic values tabulated by decay class in (Coomes et al.2002;Holdaway et al.2017) were used for all species.Finally,Carbon_Fraction was assumed to be 0.50 for all dead material.Where the species was unknown,a basic density of 477 kg·m?3was used,this being the volume-weighted mean of the tabulated density for all dead material of known species in the inventory.Carbon in stumps and dead standing spars of tree ferns,cabbage trees and palms was estimated directly from measured height and DBH by calculating the carbon in an equivalent live stem using the equation in Beets et al.(2012b) and multiplying this by the decay modifier.

Carbon stock change estimation by pool

To estimate changes in the AGB pool between two inventory measurement cycles,each measured stem was followed over time.With this approach carbon change is calculated for each stem and summed for the plot.This method requires accounting for ingrowth stems (below the DBH threshold size at Time 1) and missed DBH measurements occurring either at Time 1 or Time 2.DBHs from missing and ingrowth trees were estimated for Time 1 from measurements at Time 2 using linear regressions fitted using the SAS MIXED procedure.Models of the following form were fitted separately for each plot,predicting DBH at Time 1,Dijk1,for plot i,species j,and stem k,from DBH at Time 2 Dijk2:

where,aiand biare intercept and slope parameters,siis a random species effect assumed normally distributed,and eijkis a normally distributed error term.This approach was used for tree and shrub species.For tree ferns,palms and cabbage trees,missing DBH at Time 1 was set equal to DBH at Time 2 on the assumption that for these species,DBH does not change greatly over time(unpublished inventory data confirmed this).Although it was considered far more likely that this procedure would correct for an under-prediction of carbon at Time 1,it is acknowledged that there is some potential for the procedure to cause over-predictions of carbon at Time 1 if in fact the missing stem was actually present but misidentified with an incorrect or missing tag.The stem would then be counted twice at Time 1 leading to an underprediction of carbon stock change between the two periods.A similar procedure was used for estimating missing DBH at Time 2 from DBH at Time 1 for the small number of stems with missing DBH in the second measurement period (usually because stems were inaccessible by field crews for various reasons).

Calculation of plot-level estimates of stocks and stock changes in nested plots

The above procedures were used to estimate the carbon for each live stem and standing or fallen piece of dead wood.Estimates of carbon per hectare were then calculated for each grid location by summing the carbon estimate (kg) in each live stem and dead wood piece,divided by the horizontal area of the relevant nested plot it was sampled in,and dividing the sum by 1000 to convert to tC·ha?1.Thus,all pieces and stems ≥60 cm diameter were divided by the area of the 20 m radius subplot(0.1257 ha) while smaller stems which were only measured in the 20 m×20 m square plot were divided by the horizontal area of that plot.To ensure the strict application of the minimum diameter size for dead wood in each nested plot,minimum diameter thresholds for fallen material of 10 cm and 60 cm were applied in inner and outer plots respectively.The estimated carbon in fallen pieces lying in the inner plot with LED greater than 60 cm and SED less than 60 cm were split into two values based on an assumption of uniform taper over the length of the piece,with the carbon calculated for the portion of the length greater than 60 cm in diameter assumed to be sampled over the larger plot area,and the remainder assumed to be sampled over the inner plot area.Carbon stock change in live stems was estimated in the same way as for carbon stocks,using the sum of the change in weight of carbon in each followed live stem divided by the relevant plot area.However,no attempt was made to follow individual pieces of fallen dead material between measurements and stock change in dead wood was estimated as the difference in stocks at times 1 and 2 calculated at the plot level.

Adjustment to the above ground dead wood pool

Recent analysis by Kimberley et al.(2019) has demonstrated that the plot measurements of dead wood obtained in the inventory tend to significantly understate the true amount of dead wood in a plot.This error is thought to be due to a consistent tendency for dead wood material to be missed during the measurement process.For example,field audits conducted since 2002 revealed that dead wood stocks were being often underestimated due to unmeasurable piece shapes,loss of bark,logs and spars shattered and debris inaccessibly piled in gullies.Furthermore,under current measurement protocols,there is no attempt to measure Decay Class 4 (heavily decayed) material,nor to measure fallen stems and other dead wood material buried more than 50%in the forest floor and which is not captured anywhere else(e.g.as litter or soil carbon).There is also uncertainty around the assignment of decay classes and species identification made in the field and uncertainty about the robustness of currently used density modifiers(Kimberley et al.2019).

To estimate dead wood carbon stocks and their change we therefore followed the approach described by Kimberley et al.(2019).The dead wood measurements from the 2002–2007 inventory cycle were used as an initial starting value.The mortality between the inventory cycles was accounted for as input into the dead wood pool,while the predicted loss from decay over the period between measurements was calculated based on models developed by Garrett et al.(2019).This approach produced a higher estimate of dead wood for 2009–2014 than the actual field assessments due to the measurement error.To compensate for this,we determined an adjustment factor by an iterative procedure such that,when multiplied by the measured dead wood in both 2002–2007 and 2009–2014,the modelled estimate for 2009–2014 was the same as the adjusted measurement at that time.The adjustment factor derived using this approach was 1.808.In other words,the corrected estimates of dead wood carbon were obtained by increasing the measured values by 80.8%.

Fine debris and litter

During the first measurement cycle in 2002–2007,carbon in the litter pool was sampled and analysed in a random sample of 250 grid locations following Davis et al.(2004).This included the three sub-pools of fine woody debris (<10 cm diameter),litter and fermenting and humus material (FH).Total litter pool data (combined total of the three sub-pools) from these grid locations were used in the current study.As only a single field measure of fine debris and litter on a subset of plots was available it was not possible to estimate change in the litter pool over time.

Below ground dead wood (dead roots)

Estimates of the carbon in dead roots were not provided in earlier unpublished reports of natural forest carbon(e.g.,Beets et al.2009 or Holdaway et al.2014).While there are no direct measurements of carbon in belowground dead wood in the inventory,estimates of the below-ground dead wood carbon pool were produced based on Garrett et al.(2019).This study found roots to decay more quickly than above ground dead wood and on average the ratio of above ground to below ground decay constants was 0.76.Under the exponential decay model,this implies that if the same amount of material enters both pools at a constant rate,the total dry weight of the below-ground material will be 76% of the dry weight of the above ground material.Therefore,as the root/shoot ratio for live trees implies that the dry weight of below ground material in a tree or shrub that dies is approximately one quarter of its above ground material,the faster below ground decay rate implies that carbon in below ground dead roots will be 0.25×0.76=0.19 times the above ground dead wood carbon.Therefore,in the current study,carbon in below ground dead roots was estimated at each grid location as 0.19 times the above-ground dead wood pool.

Estimates of mean per hectare carbon stocks and stock changes

Plot measurements were obtained from 1046 grid locations,but not all grid locations were measured in both measurement cycles,and only a subset of plots were measured for litter and only in the first measurement cycle.Sample plots can be divided into five subsamples based on their measurement history which for convenience are labelled A–E (Table 1),with the combined samples measured in cycle 1 labelled F.

Table 1 Summary of measurements obtained during the first two cycles of the inventory

Estimates of per hectare carbon stocks for each pool in cycle 1 were obtained on the assumption that the sampling method adopted for this cycle was equivalent to simple random sampling of the mapped forest area.Mean stocks for AGB,BGB,and dead wood were therefore estimated using standard methods for simple random sampling (e.g.,Cochran 1977,Chapter 2),i.e.:

where,yiis per hectare carbon at Time 1 for AGB,BGB,or dead wood in plot i,and nFis the number of locations in subsample F (i.e.,nF=1036).The variance ofwas estimated using:

The estimated per hectare mean carbon stock in litter l and its variance v（l）,were estimated from li,the mean litter carbon in plot i,as in Eqs.6 and 7 except that the summation and variance was over locations measured for litter,i.e.,A U B.Total carbon for Time 1 in the AGB,BGB,dead wood,and litter pools,was estimated using:

where,sylis the sample covariance between yiand liin the combined subsamples A U B.

In the second measurement cycle,plots from the first cycle were stratified into ‘Shrub’ and ‘Forest’ strata based on species composition,and samples randomly selected from each stratum for measurement.Methods for analysing this ‘double sampling for stratification’ approach are described in Rao (1973).Estimates of carbon across the two strata in the AGB,BGB,and dead wood pools,and in the sum of these pools,were obtained using:

where,nFhis the number of samples in subsample F within stratum h,is the mean of yi,(the carbon in the relevant pool)in plot i of stratum h,and L is the number of strata (two in this case).The variance ofwas estimated using:

where,nEandare respectively the sample size and mean of yiin subsample E.The estimated variance was obtained using:

As litter was not measured at Time 2,its estimate from Time 1 was used in the estimate of total carbon in all four pools for Time 2,i.e.,

where,cov(yi,li) is the sample covariance between yiand liin subsample A.

Carbon stock changes for the AGB,BGB and dead wood pools and their sum (note that it was not possible to estimate change in the litter pool which was only measured in cycle 1) were estimated using Eq.10 withbeing the mean stock change for the relevant pool for plots in stratum h,and with variance estimated using Eq.11.

Effects of model prediction error

Estimates of carbon such as those provided in this study rely on complex allometric relationships and models used to estimate carbon from stem and log measurements,and it is desirable to account for prediction errors in these models when considering the uncertainties of carbon estimates.The variances of estimates derived from plot-to-plot sampling variation described in the previous section take account of the natural variability of tree size and distribution within the forest and allow for random measurement errors but take no account of errors in the underlying models used to estimate carbon from tree measurements.Estimates of model prediction uncertainties for carbon in each pool and for total carbon based on an analysis of the underlying models used to estimate each pool,are shown in Table 2.Details of their derivation are presented in the Supplementary Material.

Because the same models are used for predicting stocks at both Times 1 and 2,their errors are clearly not independent,and in practice may be close to identical.If these errors,expressed as percentages of each respective mean are indeed identical,then the estimate change in stocks (i.e.,stock at Time 2 minus stock at Time 1),will also have the same percentage error.For example,from Table 2,the model uncertainty for carbon stock in all pools is ±5.2% and it can be assumed that the uncertainty for stock change will therefore also be ±5.2% of itsmean.This implies,for example,that if the estimated stock change is close to zero,then its prediction error will also be close to zero.As model prediction errors are independent of errors from sampling variation,they can be combined using standard methods of combining independent variances giving the following estimates of uncertainty (see Page 3.28 in Volume 1 of Intergovernmental Panel on Climate Change (IPCC) 2006):

Table 2 Model prediction uncertainty for carbon stocks for each pool and for the sum of all pools.Uncertainties are shown as 95% confidence intervals expressed as a percentage of the mean stock

Estimates of carbon stocks and stock change by forest alliance

New Zealand’s natural forest has recently been classified into Forest Alliances based on species composition assessed in the second measurement cycle of the permanent sample plots of the natural forest inventory(Wiser 2016;Wiser and De Cáceress 2018).This classification was applied to assign the calculated carbon stock and carbon stock change to the 28 Forest Alliances identified in New Zealand (plus one unclearly classified group).The Forest Alliances were also aggregated into broad physiognomic groups following Wiser et al.(2011)with shrublands and other forests (together termed“regenerating forest”) and four forest groups (together termed“tall forests”).Note that the definition of“regenerating forest”was somewhat broader than the definition of the“Shrub”stratum used for selecting samples in cycle 2.

Results

Overall carbon stocks and stock changes

Our analysis of inventory data from permanent plots at 1046 grid locations within New Zealand’s pre-1990 natural forest shows it is in a state of carbon balance.Carbon stocks in all live and dead pools combined averaged 227.0±14.4 tC·ha?1in the first period and 227.2±14.5 tC·ha?1in the second period (Table 3).Around 73% of assessed carbon was stored in the living biomass (AGB and BGB),while dead wood and litter stored the remaining 27%.Based on its mapped area,pre-1990 natural forest in New Zealand contains 1759±112 MtC in live and dead biomass pools excluding soil carbon.The average of total carbon stocks in tall forests was approximately four times higher than carbon stocks in regenerating forest.The estimate for total carbon stocks in tall forests was 252.4 tC·ha?1for the first period and was unchanged for the second period,compared to estimates of 53.6 and 58.1 tC·ha?1in regenerating forest for the same periods.

Table 3 Estimates of carbon stocks based on 1046 plots measured in pre-1990 natural forest.Estimates are shown for all forests,and for tall forest and regenerating forest based on species composition.Plus-or-minus values are estimated 95% confidence intervals calculated using two methods;the first method is based solely on sampling variation between plots(CIs);and the second method combines the effects of sampling variation and prediction uncertainty (CIs&p)

Over the 7.7 years between the first and second measurement periods covered by the inventory,there was no significant change in total carbon across all plots,with the estimated positive carbon stock change over this period being 0.2±1.4 tC·ha?1(or 0.03±0.18 tC·ha?1·yr?1) which does not differ significantly from zero(Table 4).However,tall forests and regenerating forests which showed very different levels of carbon stocks(Table 3),also differed in carbon stock change (Table 4).In tall forest,the stock change between the two measurements was ?0.4±1.6 tC·ha?1(or ?0.06±0.20 tC·ha?1·yr?1),not differing significantly from zero,indicating that tall forests are in carbon balance.For regenerating forests carbon stock change was significantly positive,with stocks increasing by 4.9±1.9 tC·ha?1between measurements or 0.63±0.25 tC·ha?1·yr?1.Note that the difference in carbon stock estimates at Times 1 and 2 from Table 3 can differ from the stock change estimates in Table 4.This is because the stock estimates in Table 3 incorporate information from plots measured only once in either cycle whereas the estimates of stock change in Table 4 are based only on plots measured in both cycles.The stock change estimates in Table 4 are therefore the most direct and reliable estimates of stock change.Note also that because the change is close to zero,model prediction uncertainty contributions little to overall uncertainty which is dominated by uncertainty due to sampling variation.Overall,the estimates of stockchange are very precise and an annual change in carbon stocks of as little as 0.08% would have been detected as statistically significant.

Table 4 Estimates of carbon stock changes between 2002 and 2007 and 2009–2012 for all forests,and for tall forest and regenerating forest separately within the pre-1990 natural forest.Plus-or-minus values are estimated 95% confidence intervals calculated using two methods;the first method is based solely on sampling variation between plots (CIs);and the second method combines the effects of sampling variation and prediction uncertainty(CIs&p)

Carbon stocks and stock changes by forest alliance

A breakdown of carbon stocks into forest alliances based on the method described in Wiser and De Cáceres(2013) is shown in Table 5.The highest stocks were in the silver beech-red beech-kamahi forest alliance (360.5 tC·ha?1) which is within the broader class of beechbroadleaved forests (304.1 tC·ha?1).This was followed by the hard beech-kamahi forest alliance (330.6 tC·ha?1),a type of beech-forest (243.2 tC·ha?1).Forests with higher podocarp abundance such as Beech-Broadleaved-Podocarp forest and Broadleaved-Podocarp Forests had overall lower carbon stocks (260.3 and 217.6 tC·ha?1respectively).

Stock changes by alliances grouped into tall forests and regenerating forests are shown in Table 6.There were few statistically significant differences in stock changes evident among the 20 forest alliances forming the tall forest group sensu Wiser et al.(2011).The only tall forest alliance with a significant change in carbon stocks was the Kamahi-podocarp-forest which showed a significant decline in carbon stocks (?8.0±6.1 tC·ha?1).Most other tall forest alliances showed small changes in stocks that were not significantly different from zero,suggesting that tall forest alliances in general are in balance.In contrast,out of the eight alliances belonging to the regenerating forest group,three showed significant positive carbon stock changes and one,Gorse shrubland with Cabbage trees (Cordyline australis (G.Forst.)Hook.f.),dominated by the non-native Gorse (Ulex europaeus L.) showed a loss in carbon stocks between measurements.The three alliances with a significant positive carbon stock change are characterised by the strong compositional presence of Kanuka (Kunzea ericoides(A.Rich) Joy Thomps.).The remaining regenerating forest alliances did not show significant losses or gains between measurements.

Above ground live carbon:the effects of growth and mortality

Net change in above ground live carbon can be partitioned into gross increment and mortality.The net change in carbon is the sum of these two components.To calculate gross increment and mortality we utilized plots that were measured twice (n=925).The gross increment and mortality in above ground biomass were calculated overall,and separately for tall and regenerating forests (Table 7).The calculation shows that the increase in carbon in growing live stems is higher in tall forest than in regenerating forest,averaging 1.29 tC·ha?1·yr?1in tall forest and only 1.02 tC·ha?1·yr?1inregenerating forest.However,in tall forests losses in carbon from mortality are 1.33 tC·ha?1·yr?1,much higher than in regenerating forest (0.55 tC·ha?1·yr?1),more than offsetting the higher gain in carbon from tree growth.The lower level of mortality in regenerating forests (shrublands) is the reason for a statistically significant net gain in carbon in these forests despite their lower gross increment (19% less) compared to the tall forest group.It is noteworthy that the apparent small loss in carbon in tall forest between cycle 1 and 2 is not statistically significant and suggests that carbon stocks may oscillate around a long-term equilibrium level determined by the balance between growth and mortality,with the timing and intensity of storm events periodically causing departures from this.

Table 6 Estimates of carbon stock change in the AGB,BGB and dead wood pools between 2002 and 2007 and 2009–2012 in the tall forest and regenerating forest group as defined by Wiser et al.(2011) and their alliances.Changes in stock values in bold differ significantly from zero.Plus-or-minus values are estimated 95%confidence intervals calculated using two methods;the first method is based solely on sampling variation between plots (CIs);and the second method combines the effects of sampling variation and prediction uncertainty(CIs&p)

Table 6 Estimates of carbon stock change in the AGB,BGB and dead wood pools between 2002 and 2007 and 2009–2012 in the tall forest and regenerating forest group as defined by Wiser et al.(2011) and their alliances.Changes in stock values in bold differ significantly from zero.Plus-or-minus values are estimated 95%confidence intervals calculated using two methods;the first method is based solely on sampling variation between plots (CIs);and the second method combines the effects of sampling variation and prediction uncertainty(CIs&p) (Continued)

Plots with large AGB stocks are structurally complex with often large emergent trees,and these are occasionally prone to storm damage driving large mortality losses when such large trees become uprooted or break and die in storm events.Storm related natural disturbances drive the mortality related carbon stock losses in tall forests with large AGB carbon stocks.Such disturbance events create a high variability between plots in terms of the balance between increment and mortality and the measured carbon stock in AGB and dead wood (Fig.2),but on average across New Zealand’s tall forests gross increment and mortality is in balance.

Fig.2 Remeasurement plots in tall forests(n=803) showing dead wood carbon stocks(tC·ha?1)plotted against total AGB carbon stocks(tC·ha?1).Data are from the first measurement cycle 2002–2007

The AGB carbon pool change estimates confirm the expected patterns of mortality and gross increment with increasing carbon stocks across both tall and regenerating forests.Carbon losses through mortality in regenerating forest stands are generally lower than in more mature stands due to the small amount of carbon that is lost through small-tree mortality in the early differentiating phase of a forest succession.This can be shown by plotting the relationships between gross increment,mortality and net change in carbon against carbon stocks(Fig.3).Gross increment in AGB carbon increases rapidly with increasing stocks but levels off once the AGB stocks reach about 100 tC·ha?1.Above this level of carbon stocks,the annual increase in carbon produced by growing stems is remarkably constant,averaging about 1.4 tC·ha?1·yr?1.Annual losses in AGB carbon from mortality are very low when AGB stocks are low such as in early successional stages but increase steadily with increasing AGB carbon stocks.

The net annual change in AGB carbon is the sum of the gross increment and mortality.At AGB stocks below 100 tC·ha?1,growth exceeds mortality,and the net change is significantly positive.When AGB stocks are between 100 and 150 tC·ha?1,growth and mortality are in balance,and the net change does not differ significantly from zero.However,when stocks are higher than 150 tC·ha?1,mortality exceeds growth and the net change is significantly negative.This means that on average,plots with AGB carbon stocks above 150 tC·ha?1are more likely than not to lose carbon over a reasonable period of time.The relationships shown in Fig.3 explain why regenerating forest,which has a net AGB carbon stock averaging only 34.7 tC·ha?1,is gaining carbon,while tall forest,which has a net AGB carbon stock averaging 148.8 tC·ha?1is in carbon balance.

Fig.3 Regression models showing how annual increase in AGB carbon from growing stems(top dashed line with longer dashes)and annual losses in AGB carbon from mortality (lower dashed line with shorter dashes)interact to generate the annual net change in AGB carbon(solid line)in relation to mean AGB carbon stocks in New Zealand’s natural forest.For each regression model,95%confidence intervals are shown by dotted lines above and below the prediction line

Discussion

Sampling and methodological considerations

The estimates of carbon stocks and stock changes presented in this paper are the best yet provided for the vast expanse of natural forests in New Zealand.They are based on a representative forest inventory with a proven systematic sampling design with random starting location,utilise the complete plot data acquired using an appropriate permanent plot design,and include a sophisticated modelling approach to improve the prediction of dead wood.

Most national forest inventories are based on nationally representative,often random or grid-based,sampling designs covering the forested land area of the country(Tomppo et al.2010) to ensure estimates are unbiased.Plot densities are chosen to achieve a desired level of precision for key variables (Tomppo et al.2009).Toachieve full representativeness,all grid-points should be sampled during an inventory period.Or,if only a subset of plots are measured in a particular inventory period,these need to be selected at random,or using some other appropriate method (e.g.,spatially balanced sampling,Brown et al.2015)to ensure that estimates are unbiased.Especially in forest inventories of large expanses of natural forests such as in the tropics or in the boreal zones,such grid-based field inventory approaches are challenging to implement due to topography and terrain.Despite the inherent difficulties of terrain in New Zealand with natural forests found in the least accessible parts of steep mountain ranges such as the NZ Southern Alps and Fiordland,93%of the 1116 points randomly selected from the grid were measured in the first measurement cycle.During the second cycle,all previously measured points were available for remeasurement,and 912 were selected using a stratified random design,along with 10 previously unsampled points from the grid.These sampling methods would provide unbiased estimates,provided carbon stocks and stock changes in the 7% of selected locations that were not measured do not differ appreciably from the national average.Some 3% of locations could not be sampled due to inaccessibility(e.g.,due to hazardous and dangerous steep terrain),and although it is likely that carbon stocks at such locations will differ from the national average,their effect can only be minor given the small numbers of locations involved.Access was denied to a slightly larger number of points representing 4% of selected locations,but again,the effect of this omission is also likely to be minor.

Although our estimates of the total carbon stock are similar to those of the earlier study by Holdaway et al.(2017),the estimates of carbon by pool differ significantly.Table 8 shows how differences in methodology between the two studies affected the estimates of carbon stocks and stock changes by pool.The earlier study reported estimates of biomass and annual changes in biomass.For comparison,these have been converted into carbon using the same carbon fractions used in our study.

Table 7 Estimates of annual AGB increase in carbon from growth of living trees,losses from tree mortality,and net change(±95%CI)for all pre-1990 forest,and for tall forest and regenerating forest,based on plots measured in both cycles of the inventory(n=925)

Table 8 Effect of changes in ethodology on estimated carbon stocks and stock changes by pool in the current study compared with Holdaway et al.(2017)

Both studies used the allometric models of Beets et al.(2012b) to predict carbon in stems which require whole stem wood and bark density tabulated by species.Holdaway et al.(2017) used a database of wood densities tabulated for 113 species whereas we used the densities tabulated by Beets et al.(2012b).The use of different density tables was the largest source of disparity between the studies with our estimate of total carbon being decreased by 14.7 tC·ha?1compared with earlier estimate.The stem volume model requires compatible wood and bark whole stem densities which take account of the effects of bark,irregular stem shape and hollow stems.As noted by Beets et al.(2012b),these densities are lower than basic wood densities and they suggested using an adjustment of 0.905 to convert breast height outer wood density to a density compatible with the volume model.As Holdaway et al.(2017) made no such adjustment,their methods would be expected to over-predict carbon in the AGB and BGB pools.

The second largest source of disparity in the estimate of carbon stock was caused by the inclusion of data from the 20 m radius plots in our study which resulted in a reduction in estimated total carbon of 11.6 tC·ha?1.Plotdesign can play a crucial role in inventory sampling.The initial design of the New Zealand natural forest inventory (Coomes et al.2002) critically accounted for the wider spatial distribution of large trees by including a larger sample plot (20 m radius) to capture these trees accurately.Holdaway et al.(2017) used only the 20 m×20 m plots,their reason for excluding data from the 20 m radius plots being the“unreliability for New Zealand’s dense natural forests and steep,highly dissecting terrain”of such larger plots.We however followed the example of Coomes et al.(2002) who successfully tested the usability of a 20 m radius subplot in New Zealand.To add to our confidence in using the 20 m radius subplot,the continuous auditing and quality assurance programme of the natural forest inventory found no major issues in the implementation and use of a horizontal 20 m radius large tree plot over the existing three inventory cycles.Moreover,the inclusion of the 20 m radius plot data in the analysis avoids a bias in large tree stem numbers associated with using only the 20 m×20 m square plots(Paul et al.2019) which is the reason it results in a reduction in our estimate of total carbon.Our use of larger plots also aligns us with research that identified smaller plot size as an issue when monitoring forest dynamics (Clark and Clark 2000;Wagner et al.2010;McMahon et al.2019;Hetzer et al.2020).Inclusion of these additional data also provided more precise estimates of carbon stocks and stock changes.Confidence intervals of total carbon stocks and stock changes based on sampling variation were respectively 17% and 39%lower than in the earlier study.

Other methodological differences between the two studies were the adjustment to account for underestimation of the dead wood pool in our study which increased the estimate of carbon stock by 10.6 tC·ha?1,the addition of an estimate of below ground dead wood in our study which added 6.4 tC·ha?1,and differences in the estimates of the litter pool which added 8.7 tC·ha?1.The reason for the difference in litter pool estimates is unclear as both are based on the same data,but we note that our estimate is simply the mean of litter measurements from a random selection of grid points which must provide a robust and unbiased estimate.

Both studies estimated changes in carbon stocks to be close to zero and the methodological differences between the studies did not greatly affect this result.The biggest difference in stock change estimates between the studies was caused by data corrections or through differences in the height models (we were unable to separate out these two effects) with these reducing our estimate of stock change between the two measurement cycles by an average of 2.5 tC·ha?1compared to the earlier study.

Comparison with other studies

Keith et al.(2009) showed that temperate forests hold significant biomass carbon stocks and recorded the highest above ground biomass carbon levels in Australian Eucalyptus regnans forests (1819 tC·ha?1) with an averaging live above ground carbon biomass pool of 1053 tC·ha?1(67% of total above ground biomass carbon).In comparison the New Zealand forest type with the highest mean total carbon pool recorded in our inventory(first measurement cycle),silver beech-red beech-kamahi forest,averaged 207 tC·ha?1in the live aboveground carbon pool with one plot measuring 661 tC·ha?1of aboveground live biomass carbon.While we therefore did not find carbon stocks as high as those reported by Keith et al.(2009) for Australia,many of New Zealand’s forest alliances have greater above ground carbon pools than boreal and tropical forests summarised by Keith et al.(2009).This is despite the fact that our representative sampling approach did not target particular forest stands nor forest types.Based on a comparison with Goldstein et al.(2020),our average estimates of total carbon stocks for all New Zealand tall forests probably exceed the recent global estimates of temperate conifer and broadleaf forests.This comparison is difficult because we did not estimate soil carbon,while Goldstein et al.(2020)did not account for dead wood and litter pools.

Our forest classification does not rely on structural components that are correlated with carbon stocks (such as tree height or tree species stocking).The forest groups and alliances we used are only classified by species assemblage,meaning that various development stages are included in the average carbon stock for each forest type (i.e.,stands comprised of younger and smaller trees of species such as red beech regeneration following natural disturbance or windfall are included with larger trees in mature stands).We would therefore expect our carbon stock estimates to be lower than those of mature stands.

Our carbon stock estimates do compare more with northern hemisphere temperate forests such as those of the Pacific Northwest coniferous forests which average approximately 344 tC·ha?1in the above-ground tree biomass,including dead wood (Smithwick et al.2002).Our estimates are also in line with temperate South American forest types that are closely related to the New Zealand beech forest alliances (Fuscospora and Lophozonia species),such as the Nothofagus dombeyi Mirb.Oerst.and Nothofagus alpina Popp.&Endl.dominated forests of southern Chile.These forests have an average of 877 tons biomass per ha in the AGB live pool (Schlegel and Donoso 2008).

Carbon sequestration

While the importance of old growth natural forests as a reservoir of terrestrial carbon is undisputed,their potential as a carbon sink is currently under debate.A number of studies have shown the ongoing sequestration potential of large trees in old growth forests (Stephenson et al.2014),or ongoing carbon sequestration in tropical(Malhi et al.2004;Lewis et al.2009) or boreal and temperate zones (Carey et al.2001;Luyssaert et al.2008;Pan et al.2011).That old-growth forests do still sequester carbon was reasoned on an individual tree basis by Stephenson et al.(2014) as carbon accumulation increases with tree size and large trees are a feature of old growth forests.Other work however suggests that old growth forest carbon stocks decline or reach an equilibrium between growth and mortality (Xu et al.2012;Brienen et al.2015).The latter is the case for New Zealand on a national level based on our representative forest inventory data,which are not compromised by the selection of plot locations and/or focus on particular study sites (Fisher et al.2008).Although carbon stocks across the entire natural forest are approximately in balance,total carbon in regenerating forest increased by an average of 0.63 tC·ha?1·yr?1during the study period,and there is some indication of a decline in carbon stocks in the tall forest by an average of 0.05 tC·ha?1·yr?1for the same period.The estimate of change for regenerating forests is statistically highly significant,although the estimate for tall forest is not.Our estimate of the increase in carbon stocks in regenerating forest is lower than that of Holdaway et al.(2014) (1.39 tC·ha?1·yr?1),largely due to the use of the revised forest type classification system of Wiser (2016) as Holdaway et al.(2014) used an earlier classification provided by Wiser et al.(2011).If the Wiser (2016) classification is applied to the Holdaway et al.(2014) stock estimates,their estimate of annual change reduces to 1.1 tC·ha?1·yr?1.Holdaway et al.(2014) also included carbon estimated using the“Shrub”measurement method (Payton et al.2004b) which we did not include in our study,because it was found that the field assessment method for shrubs was potentially biased (Mason et al.2014).Carbon estimated using the shrub method contributed 0.18 tC·ha?1·yr?1to the change between the two measurement periods for regenerating forest in the earlier study.

We recognise that the estimation of the dead wood pool based on field measurements is difficult and therefore often imprecise in many studies.However,it is critical to understand the role of mortality on the dynamics in this pool,and also to recognise the interaction between mortality and growth in old growth forests.Figure 3 describes this interaction and shows why New Zealand’s natural forests are in carbon balance.The average AGB carbon stock across all plots in the inventory is 134 tC·ha?1,the level at which the gain in carbon from growth of living stems is precisely in balance with the loss in carbon from mortality.This suggests that there have been no large scale major external factors(e.g.national scale natural disturbances) impacting forest growth or mortality in New Zealand’s natural forests over recent decades.While disturbances do occur regionally or locally,the frequency and intensity of such effects has had a balancing effect such that re-growth more-or-less balances mortality over time nationally.

The natural forest inventory as an unbiased monitoring framework will allow further analysis to identify long-term trends in the future.This will be critical e.g.to understand the national impact of past and future management of these protected forests,specifically in regard to herbivore control or the continuous monitoring of long-term dynamics under climate change.The tagging of trees will allow the detailed analysis of mortality of specific tree species and their growth rates and ingrowth over the coming years.Moreover,when the third measurement cycle has been completed,the data quality can be expected to improve.

Conclusions

Our study presents refined methods for estimating carbon stocks and stock change in New Zealand’s natural forest,highlighting their importance as a significant terrestrial carbon pool for New Zealand and the southern hemisphere temperate biome.Based on the complete and representative forest inventory,we provide for the first time carbon stock estimates that include the full carbon stock assessment of live and dead tree biomass,including components of dead wood carbon that are very difficult to assess,achieved through recent advances in dead wood modelling.Finally,the uncertainties of estimates of stocks and stock changes provided in this paper take account of both model prediction error and uncertainty associated with sampling variation between plots,of which the latter is the dominant source of variation.

Our estimates of stock change in particular are very precise and strongly suggest that New Zealand’s natural forests are in carbon balance and are neither a carbon sink nor source,although some regenerating forests are sequestering carbon.Due to its large size the natural forest carbon pool is an important reservoir that needs to be managed to ensure that large emissions through direct or indirect human-caused forest degradation are avoided.Investment in consistent long-term monitoring will be fundamental to detect any such changes in this critical terrestrial carbon pool.

Supplementary Information

The online version contains supplementary material available at https://doi.org/10.1186/s40663-021-00312-0.

Additional file 1:Table S1.Description of height estimation method including relevant regression coefficients and estimation of prediction errors and their uncertainties.

Acknowledgements

We thank the New Zealand Ministry for the Environment for leading the LUCAS programme and all those that are involved in the design,data collection and data management.Special thanks go to the LUCAS team at the Ministry for the Environment and the dedicated field staff involved in the programme.

Authors’contributions

TP,MK and PB developed the idea and conjointly prepared the manuscript;MK designed analysis approach;TP collated the data;MK,TP and PB analysed the data;All authors contributed critically to the drafts and gave final approval for publication.

Funding

The New Zealand Ministry for the Environment provided funding to undertake data analysis and preparation of this manuscript under Statement of Work 21078.Additional support was provided by the New Zealand Ministry for Business,Innovation and Employment Core funding to Crown Research Institutes.

Availability of data and materials

The data that support the findings of this study are available from the New Zealand Ministry for the Environment,but restrictions apply to the availability of these data,which were used under license for the current study,and so are not publicly available.Data are however available from the authors upon reasonable request and with permission of the New Zealand Ministry for the Environment.

Declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Author details

1Scion,Private Bag 3020,Rotorua 3046,New Zealand.2Environmental Statistics Ltd,Whangapararoa 0930,New Zealand.3retired (previously Scion),Rotorua 3010,New Zealand.

Received:4 May 2020Accepted:20 May 2021