Modelling above-ground carbon dynamics using multi-temporal airborne lidar: insights from a Mediterranean woodland

Abstract

Woodlands represent highly significant carbon sinks globally, though could lose this function under future climatic change. Effective large-scale monitoring of these woodlands has a critical role to play in mitigating for, and adapting to, climate change. Mediterranean woodlands have low carbon densities, but represent important global carbon stocks due to their extensiveness and are particularly vulnerable because the region is predicted to become much hotter and drier over the coming century. Airborne lidar is already recognized as an excellent approach for high-fidelity carbon mapping, but few studies have used multi-temporal lidar surveys to measure carbon fluxes in forests and none have worked with Mediterranean woodlands. We use a multi-temporal (5-year interval) airborne lidar data set for a region of central Spain to estimate above-ground biomass (AGB) and carbon dynamics in typical mixed broadleaved and/or coniferous Mediterranean woodlands. Field calibration of the lidar data enabled the generation of grid-based maps of AGB for 2006 and 2011, and the resulting AGB change was estimated. There was a close agreement between the lidar-based AGB growth estimate (1.22 Mg ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$) and those derived from two independent sources: the Spanish National Forest Inventory, and a tree-ring based analysis (1.19 and 1.13 Mg ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$, respectively). We parameterised a simple simulator of forest dynamics using the lidar carbon flux measurements, and used it to explore four scenarios of fire occurrence. Under undisturbed conditions (no fire) an accelerating accumulation of biomass and carbon is evident over the next 100 years with an average carbon sequestration rate of 1.95 Mg C ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$. This rate reduces by almost a third when fire probability is increased to 0.01 (fire return rate of 100 years), as has been predicted under climate change. Our work shows the power of multi-temporal lidar surveying to map woodland carbon fluxes and provide parameters for carbon dynamics models. Space deployment of lidar instruments in the near future could open the way for rolling out wide-scale forest carbon stock monitoring to inform management and governance responses to future environmental change.

Introduction

The world's forests are currently acting as an important carbon sink, in 2000–2007 taking up 2.3 $\pm$ 0.5 PgC each year compared with anthropogenic emissions of 8.7 $\pm$ 0.8 PgC (Pan et al., 2011). For this reason, the international community recognises that forest protection could play a significant role in climate change abatement and that the feedback between climate and the terrestrial carbon cycle will be a key determinant of the dynamics of the Earth System (Purves et al., 2007). However, there is major uncertainty over forest responses to anthropogenic global change, and concerns that the world's forests may switch from being a sink to a source within the next few decades (Nabuurs et al., 2013; Ruiz-Benito et al., 2014b), through gradual effects on regeneration, growth and mortality, as well as climate-change related disturbance (Frank et al., 2015). For instance, severe droughts in many parts of the world are causing rapid change, killing trees directly through heat-stress and indirectly by fire (Allen et al., 2010). Disturbance events can cause major perturbations to regional carbon fluxes (Chambers et al., 2013; Vanderwel et al., 2013). A major goal in biogeosciences, therefore, is to improve understanding of the terrestrial vegetation carbon cycle to enable better constrained projections (Smith et al., 2013).

In this context, remote-sensing methods for modelling above-ground storage of carbon in biomass have received much recent attention, with airborne light detection and ranging (lidar) showing the most potential for accurate and large-scale applications. Lidar metrics of canopy structure are highly correlated with field-based estimates of above-ground biomass (AGB) and carbon (AGC) (Drake et al., 2003; Lefsky et al., 2002). With such relationships being repeatedly demonstrated, it has been possible to develop a conceptual and technical approach linking plot-based carbon density estimates with lidar top canopy heights using regional inputs on basal area and wood density (Asner and Mascaro, 2014). With the increasing availability of multi-temporal (repeat survey) lidar data sets, including some of national coverage, a few researchers have started to use lidar in large-scale studies of vegetation productivity and carbon dynamics (Englhart et al., 2013; Hudak et al., 2012) as well as forest disturbance and gap dynamics (Blackburn et al., 2014; Kellner and Asner, 2014; Vepakomma et al., 2008, 2010, 2011). As such, and despite its high costs, lidar is transitioning from research to practical application, notably in supporting baseline surveys and monitoring of carbon stocks required for the implementation of the REDD$+$ mechanism (Reducing Emissions from Deforestation and Forest Degradation) (Asner et al., 2013). However, monitoring carbon fluxes using multi-temporal lidar is technically challenging because instrument and flight specifications vary over time (Réjou-Méchain et al., 2015).

The applications of airborne lidar for modelling AGB and AGC have largely been tested in cool temperate and tropical forest systems (see Zolkos et al., 2013). Less attention has been given to the effectiveness of the technology for the modelling of biomass and carbon in sub-tropical and Mediterranean climate zones dominated by dry woodlands. These woodlands have lower carbon densities, but represent important global carbon stocks due to their extensiveness and also vulnerability in the face of climate change (Ruiz-Benito et al., 2014b). As elsewhere in Europe, carbon stocks in such woodlands have been increasing in recent decades (Nabuurs et al., 2003, 2010; Vayreda et al., 2012), as woodland management for charcoal and timber has declined in profitability. However, with Earth System models predicting some of the most severe warming and drying trends of anywhere in the world (Giorgi and Lionello, 2008; Valladares et al., 2014), abrupt shifts in increasing fire frequency and intensity may reverse such trends across the Mediterranean region (Pausas et al., 2008). Lidar has been used to measure carbon stocks in some Mediterranean woodlands (García et al., 2010) but, to our knowledge, not for measuring carbon dynamics.

In this study we demonstrate the potential to build a patchwork dynamics simulator for the biomass and carbon dynamics in Mediterranean woodlands based on multi-temporal lidar data (Fig. 1). Our aim is to model the direction and rate of landscape-scale AGC change for mixed oak-pine woodland in central Spain. We first calibrate a lidar top-of-canopy height model using selective ground-based estimations of tree- and plot-level biomass. The lidar-based AGB growth models are then validated using two independent data sets: the Spanish National Forest Inventory (SFI) and tree-ring measurements, before parameterising a simulation model to explore the dynamics of carbon change over a 100-year period. In doing so, we explore sensitivity of the long-term carbon sequestration potential of the regional landscape to increasing forest fire frequency, as is to be expected under future climate change.

Methodological approach.

[Figure omitted. See PDF]

Methods

Study area

Alto Tajo (40${}^{\circ }$47${}^{\prime }$ N, 2${}^{\circ }$14${}^{\prime }$ W) is a Natural Park (32 375 ha) situated in the Guadalajara province of Central Spain. The dominant woody vegetation is Mediterranean mixed woodland, comprising Pinus sylvestris, P. nigra, Quercus faginea, Q. ilex, Juniperus oxycedrus and J. thurifera. The region has a complex topography ranging from 960 to 1400 m a.s.l. The mean annual temperature here is 10.2 ${}^{\circ }$C, with mean annual rainfall of 499 mm.

Contained within the Park is one of the six Exploratory platform sites contributing to FunDivEurope: Functional Significance of Biodiversity in European Forests (Baeten et al., 2013). Field data used in the current study were taken from plots surveyed as part of this programme. The landscape-level analysis focused on a belt overlapping this area and running 20 km north–south and 3 km east–west (Fig. 2).

Study area. Shown in lighter green, mixed forest, and darker green, coniferous forest. Other land covers (including agricultural) in shades of grey, with darkest grey indicating an area burned by forest fire in 2005 and excluded from these analyses. The three north-south parallel strips show the lidar survey coverage.

[Figure omitted. See PDF]

Plot-based tree measurements and allometric biomass modelling

Field measurement of plots was undertaken in March 2012. Each plot was of dimension 30 $\times$ 30 m and was carefully geo-located, recording GPS corner coordinates and orientation using a Trimble GeoXT – Geoexplorer 2008. Measurements were made of trees and shrubs of diameter at breast height (DBH) > 7.5 cm, given that smaller sizes contribute less to plot-level biomass (Stephenson et al., 2014). The following were measured and recorded: position within plot, species, height, height of lowest branch, DBH (at 1.3 m), and crown diameter (two orthogonal measurements). A vertex hypsometer was used for the crown dimensions.

The above-ground biomass of individual trees was estimated according to published allometries, and summed to arrive at plot and hectare totals. The allometric equations of Ruiz-Peinado et al. (2011) and Ruiz-Peinado et al. (2012) were used for softwood species (Juniperus and Pinus) and hardwood species (Quercus), respectively (Appendix A). The equations were developed from tree samples across Spain including sites close to the Alto Tajo study area. The equations for Juniperus thurifera were applied to the other two junipers (J. oxycedrus and J. phoenicia) as well as box (Buxus sempervirens). In all cases, the equations compartmented the biomass into trunks and large, medium and fine branches and/or leaves, using DBH and tree height data.

Lidar surveys, calibration and above-ground biomass and carbon change analysis

The lidar surveys were undertaken by the NERC Airborne Research and Survey Facility (ARSF) and took place on 16 May 2006 (project WM06_04; García et al., 2011, 2010) and 21 May 2011 (project CAM11_03). A Dornier 228 aircraft was employed for both, but lidar instruments differed between years: Optech ALTM-3033 in 2006 and Leica ALS050 in 2011. Instrument and flight parameters are given in Table 1. Simultaneous GPS measurement was carried out on the ground allowing for differential correction during post-processing.

Specifications for the lidar surveys undertaken at Alto Tajo (Spain) in 2006 and 2011.

	2006	2011
Lidar sensor	Optech-ALTM3033	Leica ALS050
Wavelength (nm)	1064	1064
Beam divergence (mrad)	0.20	0.22
Vertical discrimination (m)	4.6	2.8
Detection system	Two return	Four return
Date of deployment	16 May 2006	21 May 2011
Pulse rate frequency (MHz)	33.33	67.2–74.4
FoV (degrees)	12	40
Scan frequency (Hz)	42.4	35.8–40.0
Point density (m${}^{-\mathrm{2}}$)	0.5	2
Number of flight lines	3(N–W)	4 (E–W) $+$ 3(N–W)
Altitude (m a.s.l.)	2063–2073	2097–2140

We assumed accurate georeferencing of the 2006 and 2011 data sets during post-processing, and did no further co-registration. We performed initial modelling of terrain and canopy heights from the 2006 and 2011 lidar data sets using “Tiffs” 8.0: Toolbox for Lidar Data Filtering and Forest Studies, which employs a computationally efficient, grid-based morphological filtering method described by Chen et al. (2007). Outputs included filtered ground and object points, as well as digital terrain models (DTM) and canopy height models (CHM). The subsequent GIS and statistical analyses described below were undertaken in ArcInfo 10.0 (undertake in ArcInfo 10.0 by ESRI) and R 2.13.1 (R Development Core Team, 2011), respectively.

Spatially overlaying the lidar data set with land cover information derived from the 2006 CORINE map (EEA, 1995), indicated the local presence of two main forest types: coniferous and mixed (oak-juniper-pine) woodland. For the purposes of calibrating the lidar height models based on field-estimated biomass, only the latter forest type was adequately sampled (13 plots), so subsequent analysis and modelling focused on these mixed woodland systems. We predicted biomass as a function of top-of-canopy heights, which has been found to be a good predictor (Asner et al., 2013). Digitised plot boundaries for the 13 FunDiv plots of square 30 $\times$ 30 m were used to extract mean top-of-canopy height values from the lidar CHM (TCH${}_{\mathrm{L}}$). Reassuringly, these values were remarkably similar to the mean canopy height estimated from plot data (TCH${}_{\mathrm{P}}$), calculated from height and crown area of each tree obtained by allometric formulae (see Kent et al. 2015); there was almost a $\mathrm{1}:\mathrm{1}$ relationship between the two estimates of height: TCH${}_{\mathrm{G}}=1.79+0.999$ $\times$ TCH${}_{\mathrm{L}}$ ($R^{\mathrm{2}}=0.88$). Field-estimated AGB was modelled on the basis of lidar mean height by linear regression of log transformed variables. Our selected model (log(AGB) $=$ 3.02 $+$ 0.89 $\times$ log(TCH${}_{\mathrm{L}})$, $R^{\mathrm{2}}=0.53$, RMSE $=$ 0.28) was back-transformed and multiplied by a correction factor (CF) to account for the back-transformation of the regression error (Baskerville, 1972); the correction factor is given by CF $=$ $e^{\text{MSE/2}}$, where MSE is the mean square error of the regression model.

We used the regression model and lidar data set to map biomass and biomass change. We aggregated canopy heights at 1 m resolution to mean values per 30 $\times$ 30 m grid cell, to reduce mismatches with the field inventory plots (Réjou-Méchain et al., 2015). The aggregation was also effective in dealing with gappiness noted in the 2006 data set due to uneven distribution of scan lines and lower point density (Table 1). Negative values caused by occasional inaccuracies evident in the DTM models, especially for 2006, were removed from the data set to avoid anomalies. For each grid cell along the three north–south transects, we were able to use the mean height–AGB regression relationship to generate estimates of AGB in 2006 and 2011, and AGB change 2006–2011.

Validation

Due to the relatively low number of ground truth plots, it was especially important to validate the lidar-modelled AGB estimates, and this was done using two different data sets. Firstly, equivalent estimates of AGB and AGB change were developed using detailed tree measurements from the Spanish National Forest Inventory (SFI). The SFI covers the forested areas of the country on a 1 km${}^{\mathrm{2}}$ grid (Villanueva, 2004). A subset of 234 SFI plots surrounding the study area and of comparable topography and climate were selected, and the data extracted for the second and third surveys (2SFI, 1992–94 and 3SFI, 2003–2006; i.e. an 11-year interval for this region). For each plot, plot-level AGB was calculated by applying the allometric equations of Ruiz-Peinado et al. (2011, 2012; Appendix A) to individual tree height and stem diameter measurements and summing these up to the plot level. Information on topoclimate (altitude, rainfall, temperature; Gonzalo, 2008) and management and/or fire disturbance were also available per plot, although areas significantly burned after the first inventory were removed from the data set.

Secondly, plot-level above-ground wood productivity values were calculated from tree-ring measurements from the same FunDiv plots used to calibrate the lidar data, according to a four-step procedure described in Jucker et al. (2014): measuring growth increments from wood cores, converting diameter increments into biomass growth, modelling individual tree biomass growth, and scaling up to plot level. For the coring, bark-to-pith increment cores were collected for a subset of trees in each plot (using a 5.15 mm diameter increment borer, Haglöf AB, Sweden). Following a size-stratified random sampling approach, one core was extracted from each selected tree at a height of 1.3 m off the ground; 12 trees per plot were cored in monocultures and 6 trees per species were cored in mixtures (Jucker et al., 2014). In this approach, plot level estimates were based on the growth of trees present in 2011 and did not account for the growth of trees that died between 1992 and 2011.

Biomass growth estimation and simulation modelling

Plotting the 30 $\times$ 30 m pixel-level AGB estimates from 2006 versus 2011 revealed a small number of outliers of AGB change that may have resulted from anomalies in the DTM and top-of-canopy modelling (see discussion). We used robust regression to remove these outliers in order to obtain reliable estimates of mean growth and its uncertainty. This was performed with the rlm command in the MASS package of R, which uses iterative re-weighted least squares (M-estimation) (Venables and Ripley, 2002). Robust regression assigns lower weights to outliers than to points close to the regression line (in our case, using a bisquare weighting function), and then uses these weights to downplay the importance of these outliers in the linear regression. On inspection of the weights, we observed that all the obvious outliers had been assigned a weight of zero, so were easily filtered out. Some 3.3 % of the data were trimmed in this way. The residuals of the remaining data set were close to normally distributed. Change in AGB was calculated for each plot in the trimmed data set as (AGB${}_{2011}$– AGB${}_{2006})/\mathrm{5}$, and the mean and standard deviation estimated. There was significant spatial auto-correlation of AGB${}_{2006}$ values (Moran's I $=$ 0.138, $p\mathit{<}\mathrm{\hspace{0.17em}0.001}$) and also AGB change (Moran's I $=$ 0.038, $p\mathit{<}\mathrm{\hspace{0.17em}0.001}$). However, following the conclusion of Hawkins et al. (2007) that regression estimates are not significantly affected by spatial autocorrelation, we considered it unnecessary to subsample the gridded data set to avoid it.

The trimmed data set was used to model AGB growth as a function of biomass, using Bayesian inference, and to create a woodland dynamics simulator. The growth model was $$\begin{array}{ll}{\text{AGB}}_{2011}=& a+b\times {\text{AGB}}_{2006}+\mathit{\epsilon }\text{where}\mathit{\epsilon }\\ \sim & \text{N}(\mathrm{0},c+d\times {\text{AGB}}_{2006})\end{array}$$ where $a$, $b$, $c$ and $d$ are parameters calculated using STAN (STAN Development Team, 2014), a Bayesian inference package. We used uninformative prior and a burn-in of 5000 iterations (well in excess of that needed for convergence), then took 100 samples from the posterior distribution. We also fitted a model containing a quadratic biomass term, but the 95 % confidence intervals of the quadratic term overlapped with zero, indicating no support for its inclusion.

Parameter values drawn from the posterior distribution were fed into a simple simulation model. We created a 5000 cell “landscape” with starting biomass sampled randomly from AGB${}_{2006}$. For each cell the annual biomass increments were estimated by drawing parameters randomly from the posterior distribution $$\mathrm{\Delta }\text{AGB}=(a+(b-\mathrm{1})\times \text{AGB}+\mathit{\epsilon })/\mathrm{5},$$ where $\mathit{\epsilon }$ was drawn at random from $N$ (0, $c+d\times$ AGB). The biomass of each cell was then altered by $\mathrm{\Delta }$AGB and the iterative process continued for 100 years. Mean AGB values for the landscape each year were recorded and plotted with 95 % confidence intervals.

We also included the effect of various fire scenarios on mean biomass change and carbon dynamics in a simplistic way. We assumed that the probability of a cell being destroyed by fire, $p$, did not depend on that cell's AGB and did not vary among years. For each time step and pixel, we decided whether a fire event had occurred in a cell by drawing random numbers from the binomial distribution, with the AGB being reset to zero as a result of a fire event. An annual probability of fire occurrence for the region of Guadalajara, based on areas burned each year from 1991 to 2010 (Ministerio de Agricultura, 2002, 2012) is $p=0.002$, whilst that from a model parameterized from topoclimatic data from southern Spain is $p=0.004$ (Purves et al., 2007). A five-fold increase in area burned as a result of a high emission climate scenario is predicted for similar forest types in Portugal (see Carvalho et al., 2009). Thus, as well as the no-fire scenario, we tested the three fire probabilities of $p=0.002$, 0.004 and 0.01 to look at the sensitivity of carbon accumulation in the mixed woodlands to a realistic range of fire frequencies. Carbon sequestration potential (mean carbon storage in biomass over the simulation period, Mg ha${}^{-\mathrm{1}}$) was calculated using the IPCC default 0.47 carbon fraction (McGroddy et al., 2004), and scaled up to a total value of carbon (and CO${}_{\mathrm{2}}$ equivalent, 3.67 $\times$ C, Mt) for all mixed woodland in the autonomous community of Castilla La Mancha (181 000 ha) under the no-fire and three fire scenarios. We acknowledge that the simulation model is basic, and since it is not spatially explicit it makes no consideration of landscape connectivity. However, the results provide insight into the likely effect of varying fire rates on carbon dynamics.

Results

Lidar estimated mean AGB of mixed woodlands was 41.8 in 2006 and 47.9 Mg ha${}^{-\mathrm{1}}$ in 2011. Mean biomass change in this 5-year period was 1.22 Mg ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$, with a considerable degree of variation around this estimate (SD $=$ 1.92 Mg ha${}^{-\mathrm{1}}$) and a large number of pixels losing biomass (Fig. 3), presumably as a result of disturbance. There was very good agreement between above-ground biomass estimated from the lidar modelling and Spanish National Inventory plots for mixed oak-juniper-pine woodland (Table 2). The lidar-based estimate is also in reasonable agreement with that calculated from the 2006 data set in an earlier analysis: 44.7 Mg ha${}^{-\mathrm{1}}$ for holm oak woodland (García et al., 2010). AGB change as modelled by the lidar approach was also close to estimates derived from the SFI and the Fundiv tree ring data (Table 2). The standard deviation of the lidar-based AGB change estimate is relatively high, probably as a result of lidar sampling and/or processing errors that are greater than measurement errors associated with plots and tree rings. From the lidar data set, there was a statistically significant but minor effect on AGB change of altitude (range 908–1322 m; $\mathrm{\Delta }$AGB $=$ 21.17–0.01 $\times$ altitude, $R^{\mathrm{2}}=0.0180$, $p\mathit{<}0.001$) and aspect (calculated as folded aspect $|$aspect–180$|$; $\mathrm{\Delta }$AGB $=$ 3.31–0.03 $\times$ aspect, $R^{\mathrm{2}}=$ 0.0057, $p\mathit{<}0.001$).

Scatterplot of above-ground biomass (AGB) estimates for 2006 and 2011: lidar (black dots), Spanish Forest Inventory (red bordered circles), with one-to-one line (black) and fitted model (green).

[Figure omitted. See PDF]

Simulation model results for AGB over a 100-year period without fire (a) and at annual fire probability of occurrence of $p=0.002$ (b), 0.004 (c) and 0.01 (d). Figures show mean (black line) and 95 % confidence intervals (grey shading).

[Figure omitted. See PDF]

Comparison of the lidar modelling of above-ground biomass (AGB) and biomass change (AGB change) with forest inventory and tree-ring data: values given are mean (and standard deviation in parentheses).

	Lidar data	Forest inventory	Tree-ring data
		data
AGB (Mg ha${}^{-\mathrm{1}}$)	41.80 ($\pm$25.68)	42.8 ($\pm$52.7)	–
AGB change (Mg ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$)	1.22 ($\pm$1.92)	1.19 ($\pm$1.17)	1.13($\pm$0. 54)
Sample size	9136 grid cells	66 plots	13 plots

Biomass change was modelled according to the relationship: $$\begin{array}{ll}{\text{AGB}}_{2011}=& 3.98+1.05\times {\text{AGB}}_{2006}+\mathit{\epsilon }\text{where}\mathit{\epsilon }\\ & \sim \text{N}(\mathrm{0},4.32+1.10\times {\text{AGB}}_{2006}.\end{array}$$

With $b=1.05$ (i.e. > 1), the woodlands are accumulating biomass over time, though the variance term is large and so some cells are losing biomass (Fig. 3). The disturbance-free simulation model showed a strong increase in accumulated AGB over the whole 100-year period (Fig. 4a). The mean AGB rose from 42.6 ($\pm$5.6) to 236.9 ($\pm$18.5) Mg ha${}^{-\mathrm{1}}$, which equates to a mean carbon flux of 1.95 MgC ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$. By modelling the occurrence of fire at probabilities of $p$ $=$ 0.002, 0.004 and 0.01, we showed its potential impact on biomass and therefore carbon accumulation (Fig. 4, Table 3). Mean (and standard deviation) values for AGB after 100 years were 200.6 ($\pm$21.1), 174.2 ($\pm$22.7), and 114.1 ($\pm$21.5) Mg ha${}^{-\mathrm{1}}$ for a fire probability of 0.002, 0.004 and 0.01 (or return rate of 500, 250 and 100 years), respectively. The effects of increasing fire occurrence also have dramatic effects on the carbon sequestration potential of the mixed woodlands considered at a regional level (i.e. Castilla la Mancha, Table 3), with the most severe fire regime reducing that potential by almost a half.

Discussion

Here we provide a demonstration of the potential of lidar remote sensing to deliver large-scale high-fidelity maps of above-ground biomass and carbon dynamics. Our lidar-based biomass growth model, estimating a mean annual growth of 1.22 MgC ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$, is in excellent agreement with the estimate independently derived from the Spanish National Forest Inventory (1.19 MgC ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$). Even though there is a large standard deviation around our estimate, the enormous sample size (9136 pixels) means that standard errors become miniscule, so our landscape level projections are delivered with high precision and reliability (Coomes et al., 2002). The number of field sampling plots used to calibrate the lidar top-of-canopy model is statistically enough given the parameters calculated and, therefore, for the purposes of our study. The coefficient of determination of the resulting model ($R^{\mathrm{2}}=0.53$) can be compared with a value of 0.67 obtained by García et al. (2010) for the same region. The difference could be due to that fact that García et al. (2010) included more plots across a greater range of woodland types, heights and carbon densities.

Average above-ground biomass (AGB) and carbon sequestration potential over a 100-year period for the four forest fire scenarios (no fire and at annual fire probability of occurrence of $p=0.002$, 0.004 and 0.01), scaled up to the regional level (181 000 ha of mixed forest in Castilla la Mancha) for carbon and carbon-dioxide equivalence.

Fire	AGB	Carbon	Regional	Regional CO${}_{\mathrm{2}}$
scenario	(Mg ha${}^{-\mathrm{1}}$)	sequestration	carbon	equivalent
		potential	(Kt)	(Kt)
		(Mg ha${}^{-\mathrm{1}}$)
No fire	124.9	58.7	10.6	39.0
$P=0.002$	111.6	52.4	9.5	34.8
$P=0.004$	101.9	47.9	8.7	31.8
$P=0.01$	77.7	36.5	6.6	24.3

In the Anthropocene era of rapid climate and environmental change, there is an urgent need for reliable large-scale monitoring of above-ground biomass and carbon stocks in forests and woodlands (Henry et al., 2015), and developing our understanding of how carbon stocks will change in the future. Forests serve the critical function of sequestering atmospheric carbon and reducing the potential rate of climate change. However, they also provide other highly important services, including provision of timber, food and other non-timber products, regulation of water cycle and habitat for biodiversity (Gamfeldt et al., 2013; Ojea et al., 2012; WRI, 2005). The amount of biomass in forest is a metric relevant to all of these functions, with an especially close relationship with sequestered and stored carbon (Boisvenue and Running, 2006). In the context of climate change mitigation and emissions target agreements made at national level, robust methodologies are needed for the regular assessment of carbon stocks in forests (Gibbs et al., 2007).

Our work demonstrates one such robust approach that has delivered a credible model of landscape-level carbon stocks and fluxes based on a 5-year interval repeat-survey lidar data set. The methodology involved identifying and discarding a small number of outliers in the AGB estimates, and it is worth reflecting on their origin. One of the challenges of multi-temporal lidar analyses are when different instruments and specifications are used in the surveys. In our case, the 2006 lidar survey had a much lower point density than for 2011, and inspection of the resulting point cloud indicated a considerably uneven distribution of the scan lines. The accuracy of the resulting terrain and canopy models will therefore be lower, potentially giving rise to some of the anomalies in our results. We sought to quantify the source of this error by performing a comparison of top-of-canopy height (TCH) models from crossing flight-lines (data not given) for both years at the 30 m grid scale, for which the standard deviation for 2006 was more than double that for 2011. TCH is known to be quite robust across different instruments (Asner and Mascaro, 2014), being less susceptible to differences in laser canopy penetration than mean canopy height (MCH) (Næsset, 2009). We considered that the size of our plots was sufficient for calibrating the system, though in comparison with larger plots: (1) errors caused by spatial misalignment of plots and lidar data are greater (Asner et al., 2009); (2) integrating measurements provides a less representative average (Zolkos et al., 2013); and (3) disagreement in protocol between lidar and field observations is greater (influenced by the effects of bisecting tree crowns in lidar data versus calling a tree “in” or “out” of the plot in field data; Mascaro et al., 2011). With regard to the latter issue, the potential error is affected by the average crown size relative to plot dimensions, such that it will be less in our situation (as it also is for boreal forest, Næsset et al., 2011), than it would be for tropical forests.

At the extensive spatial scales required, remote-sensing methodologies offer the only practicable approach to the challenge of forest monitoring, with lidar being the remote-sensing instrument of choice given its potential to characterise the three dimensional structure of canopies and understories to a high degree of accuracy and resolution. Whilst spatial and temporal lidar coverage of the terrestrial and wooded surface of the planet is still limited, and the costs still high, this situation is improving continuously. A number of national surveys have been undertaken or commissioned, and building on the experience of the GLAS (Geoscience Laser Altimetry System) instrument on ICESAT (2003–2010), the GEDI Lidar space-borne facility is planned for deployment in 2019 (Dubayah et al., 2014). With these advancements, it is an important time to develop proof of principle of lidar monitoring of forest biomass and carbon stocks and fluxes. In this respect, a number of important multi-temporal lidar studies have emerged. Typical of these are an analysis of AGB dynamics, tree growth and peat subsidence in peat swamp forests of Central Kalimantan, Indonesia 2007–2011 (Boehm et al., 2013; Englhart et al., 2013), biomass changes in conifer forests of northern Idaho 2003–2009 at the pixel, plot and landscape level and looking at the impacts of logging (Hudak et al., 2012), studies of canopy gap dynamics (Blackburn et al., 2014; Vepakomma et al., 2008, 2010, 2011), and treefall rates and spatial patterns in a savanna landscape 2008–2010 (Levick and Asner, 2013). A study employing four lidar surveys between 2000–2005 established an optimum interval (3 years) for measuring tree growth in red pine forests at an acceptable level of uncertainty (Hopkinson et al., 2008).

Our study makes an important additional contribution to this literature. It demonstrates how sampling a woodland system with a small number of field plots can effectively calibrate a lidar data set to scale up credible estimates of AGB and AGC at the landscape level. It is also novel in studying these dynamics within a Mediterranean environment. Much focus of lidar-based biomass modelling has been on tropical forest systems, given their importance to the global carbon cycle. Mediterranean woodlands hold a much lower carbon density, yet are valuable carbon stores given their extensive nature not just in the Mediterranean Basin but also other similar climate regions in the world. Furthermore, the potential effects of climate change in Mediterranean woodlands are suggested to be particularly strong (Benito-Garzón et al., 2013; Ruiz-Benito et al., 2014b). In the absence of fire in one such region, our simulation suggests a significant AGB increase from 42.6 to 236.9 Mg ha${}^{-\mathrm{1}}$ over a 100-year period (equivalent to 1.94 MgC ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$). Pan et al. (2011) estimates an annual increase of 1.68 MgC ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$ in European temperate forests in 2000–2007, whilst the annual carbon sink in Mediterranean pine plantations range between 1.06–2.99 MgC ha${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$ depending on species and silvicultural treatment (Bravo et al., 2008). Estimates provided by Ruiz-Benito et al. (2014) range from 0.55 (sclerophyllous vegetation) to 0.73 (natural pine forest) and 1.45 (pine plantation). Our own estimate of carbon sequestration potential equates to a regional carbon sequestration potential of over 10 M kg (19 kt CO${}_{\mathrm{2}}$ equivalent) for mixed woodlands in Castilla la Mancha. Such a figure can be set in the context of national level commitments to the reduction of greenhouse gas emissions of 10 % against the Kyoto base year value of 289.8 Mt CO${}_{\mathrm{2}}$ equivalent (EEA, 2014). Under Spain's `Socioeconomic Plan of Forest Activation', land use, land use change and forestry (LULUCF) is projected to absorb 20–30 Mt CO${}_{\mathrm{2}}$ equivalent per year.

The contribution of Mediterranean forests to the greenhouse gas balance sheet is vulnerable to the effects of climate change, for which the Mediterranean is a hotspot region (Giorgi and Lionello, 2008; Lindner et al., 2010). One of the mediating drivers is forest fire risk. We found that an increase in fire probability from 0.002 to 0.01 (return rate increase from 500 to 100 years) dramatically altered the carbon sequestration potential of the landscape, with carbon stocks much reduced after 100 years with the highest fire probability scenario. It is worth noting in this respect that our modelled range of fire probabilities are conservative compared to estimates used in other simulations for similar regions (e.g. 0.01–0.2 for Catalonia, Lloret et al., 2003). However, it is also necessary to note that our simplistic modelling of fire, using a set probability of a burn irrespective of factors such as landscape position and temporal variability, means that our results can only be treated as indicative of the scale of effect of different scenarios on the landscape carbon dynamics. For example, our modelling does not account for the way in which small changes in temperature and rainfall regimes could lead to tipping points of much higher risk and frequency, if not severity, of burns (Moritz et al., 2012), and dramatically different carbon dynamics outcomes.

Our modelling is neither able to account for ecophysiological factors. Tree physiology is responsive to changing temperature and soil water availability, influencing rates of regeneration, growth and mortality (Choat and Way, 2013; Choat et al., 2012; Frank et al., 2015; Williams et al., 2012). One study of low productivity forests (including Alto Tajo as a continental Mediterranean study area) showed how leaf respiration rates, and their ability to acclimate to seasonal changes in the environment, have a profound effect on whether trees can maintain productivity – and continue to act as carbon sinks – in dryland areas (Zaragoza-Castells et al., 2008).

Nevertheless, our modelling approach shows considerable promise for understanding the effects of different drivers on vegetation dynamics and making informative future predictions (Chambers et al., 2013; Coomes and Allen, 2007; Espírito-Santo et al., 2014). We compared no-fire with three different fire scenarios, but it would be equally possible to develop our approach further to consider other environmental and ecological drivers of the AGB and AGC dynamics, including tree diversity (Jucker et al., 2014; Ruiz-Benito et al., 2014a) and competition effects (Ruiz-Benito et al., 2014a, b; Vayreda et al., 2012). With regard to understanding the landscape-level carbon dynamics of Spanish forests, in further work we propose coverage of a full range of different forest types and the development of more sophisticated climate change scenarios using models based on meteorological data, environmental parameters and different IPCC projections. More widely, the further development and testing of these methods is critical for exploring the prospects for, and contribution of, forests in the global carbon cycle under future environmental change.

Allometric equations used in the estimation of tree biomass from height and stem diameter measurements (Ruiz-Peinado et al., 2011, 2012).

Pinus nigra Arn. $$\begin{array}{cccc}& \text{Stem}{W}_{\mathrm{s}}=& & 0.0403\cdot d1.838\cdot h^{0.945}\\ & \text{Thick branches}{W}_{\text{b2–7}}=& & \text{If}d\le 32.5\mathrm{cm}\text{then}Z=\mathrm{0};\\ & & & \text{If}d>32.5\mathrm{cm}\text{then}Z=\mathrm{1};\\ & & & [0.228\cdot (d-32.5{)}^{\mathrm{2}}]\cdot Z\\ & \text{Medium branches}{W}_{\text{b2–7}}=& & 0.0521\cdot d^{\mathrm{2}}\\ & \text{Thin branches}+\text{needles}{W}_{\mathrm{b}\mathrm{2}+\mathrm{n}}=& & 0.0720\cdot d^{\mathrm{2}}\\ & \text{Roots}{W}_{\mathrm{r}}=& & 0.0189\cdot d^{2.445}\end{array}$$

Pinus sylvestris L. $$\begin{array}{cccc}& \text{Stem}{W}_{\mathrm{s}}=& & 0.0154\cdot d^{\mathrm{2}}\cdot h\\ & \text{Thick branches}{W}_{\text{b2–7}}=& & \text{If}d\le 37.5\mathrm{cm}\text{then}Z=\mathrm{0};\\ & & & \text{If}d>37.5\mathrm{cm}\text{then}Z=\mathrm{1};\\ & & & [0.540\cdot (d-37.5{)}^{\mathrm{2}}-0.0119\\ & & & \cdot (d-37.5{)}^{\mathrm{2}}\cdot h]\cdot Z\\ & \text{Medium branches}{W}_{\text{b2–7}}=& & 0.0295\cdot d^{2.742}\cdot h^{-0.899}\\ & \text{Thin branches}+\text{needles}{W}_{\mathrm{b}\mathrm{2}+\mathrm{n}}=& & 0.530\cdot d^{2.199}\cdot h^{-1.153}\\ & \text{Roots}{W}_{\mathrm{r}}=& & 0.130\cdot d^{\mathrm{2}}\end{array}$$

Juniperus thurifera L. (applied for all Juniperus) $$\begin{array}{cccc}& \text{Stem}{W}_{\mathrm{s}}=& & 0.0132\cdot d^{\mathrm{2}}\cdot h+0.217\cdot d\cdot h\\ & \text{Thick branches}{W}_{\text{b2–7}}=& & \text{If}d\le 22.5\mathrm{cm}\text{then}Z=\mathrm{0};\\ & & & \text{If}d>22.5\mathrm{cm}\text{then}Z=\mathrm{1};\\ & & & [0.107\cdot (d-22.5{)}^{\mathrm{2}}]\cdot Z\\ & \text{Medium branches}{W}_{\text{b2–7}}=& & 0.00792\cdot d^{\mathrm{2}}\cdot h\\ & \text{Thin branches}+\text{needles}{W}_{\mathrm{b}\mathrm{2}+\mathrm{n}}=& & 0.273\cdot d\cdot h\\ & \text{Roots}{W}_{\mathrm{r}}=& & 0.0767\cdot d^{\mathrm{2}}\end{array}$$

Quercus faginea $$\begin{array}{cccc}& \text{Stem}{W}_{\mathrm{s}}=& & 0.154\cdot d^{\mathrm{2}}\\ & \text{Thick branches}{W}_{\text{b2–7}}=& & 0.0861\cdot d^{\mathrm{2}}\\ & \text{Medium branches}{W}_{\text{b2–7}}=& & 0.127\cdot d^{\mathrm{2}}-0.00598\cdot d^{\mathrm{2}}\cdot h\\ & \text{Thin branches}+\text{leaves}{W}_{\mathrm{b}\mathrm{2}+\mathrm{l}}=& & 0.0726\cdot d^{\mathrm{2}}-0.00275\cdot d^{\mathrm{2}}\cdot h\\ & \text{Roots}{W}_{\mathrm{r}}=& & 0.169\cdot d^{\mathrm{2}}\end{array}$$

Quercus ilex $$\begin{array}{cccc}& \text{Stem}{W}_{\mathrm{s}}=& & 0.143\cdot d^{\mathrm{2}}\\ & \text{Thick branches}{W}_{\text{b2–7}}=& & \text{If}d\le 12.5\mathrm{cm}\text{then}Z=\mathrm{0};\\ & & & \text{If}d>12.5\mathrm{cm}\text{then}Z=\mathrm{1};\\ & & & [0.0684\cdot (d-12.5)\mathrm{2}\cdot h]\cdot Z\\ & \text{Medium branches}{W}_{\text{b2–7}}=& & 0.0898\cdot d^{\mathrm{2}}\\ & \text{Thin branches}+\text{leaves}{W}_{\mathrm{b}\mathrm{2}+\mathrm{l}}=& & 0.0824\cdot d^{\mathrm{2}}\\ & \text{Roots}{W}_{\mathrm{r}}=& & 0.254\cdot d^{\mathrm{2}}\end{array}$$ Notes:$W_{\mathrm{s}}$: Biomass weight of the stem fraction (kg);$W_{\mathrm{b}\mathrm{7}}$: Biomass weight of the thick branches fraction (diameter larger than 7 $\mathrm{cm}$) (kg); $W_{\text{b2–7}}$: Biomass weight of medium branches fraction (diameter between 2 and 7 $\mathrm{cm}$) (kg); $W_{\text{b2}+\mathrm{l}}$: Biomass weight of thin branches fraction (diameter smaller than 2 $\mathrm{cm}$) with leaves (kg); $W_{\mathrm{r}}$: Biomass weight of the belowground fraction (kg); $d$: diameter at breast height (cm); $h$: tree height (m).

The project was conceived by D. Coomes and W. Simonson. Lidar analysis and first manuscript drafting was undertaken by W. Simonson. D. Coomes designed the statistical approach, and P. Ruiz-Benito provided the independent validation data and analysis based on the Spanish National Forest Inventory. F. Valladares oversaw field data collection, and with all authors contributing to the finalisation of the manuscript.

Acknowledgements

Field data were collected by T. Jucker and partners from University Stefan cel Mare of Suceava (Romania) and National Museum of Natural Sciences, Madrid (Spain). Biomass estimates were calculated by T. Jucker. The authors would like to acknowledge the personnel of the Airborne Research and Survey Facility (NERC). We thank the MAGRAMA for granting access to the Spanish Forest Inventory. WS was funded by FunDivEurope and the Isaac Newton Trust. PRB was supported by The International Post doc Fellowship Programme in Plant Sciences (PLANT FELLOWS). Edited by: U. Seibt