Use of Aerial Laser Scanning to Assess the Effect on C Sequestration of Oak ( Quercus ilex L. subsp. ballota [Desf.]Samp ‐ Q. suber L.) Afforestation on Agricultural Land

: Conversion of agricultural lands to forest plantations to mitigate rising atmospheric carbon dioxide (CO 2 ) has been proposed, but it depends on accurate estimation of the on ‐ site carbon (C) stocks distribution. The use of aerial laser scanning (ALS) data is a rapidly evolving technology for the quantification of C stocks. We evaluated the use of allometric models together with high ‐ density ALS data for the quantification of biomass and soil C stocks in a 14 ‐ year ‐ old Quercus ilex and Q. suber plantation in Southwestern Spain. In 2010, a field survey was performed and tree dasometric and biomass variables were measured. Forty ‐ five soil profiles (N = 180 soil samples) were taken systematically and the soil organic C content (SOC) was determined. Biomass and soil organic C values were regressed against individual dasometric variables and total tree height was used as a predictor variable. Aerial laser scanning data were acquired with a point density of 12 points m − 2 . Relationships among ALS metrics and tree height were determined using stepwise regression models and used in the allometric models to estimate biomass and SOC C stocks. Finally, a C stock map of the holm ‐ cork oak cover in the study area was generated. We found a tree total biomass of 27.9 kg tree − 1 for holm oak and 41.1 kg tree − 1 for cork oak. In the holm oak plantation, the SOC content was 36.90 Mg ha − 1 for the layer 0–40 cm (SOC 40 ) under the tree crown and 29.26 Mg ha − 1 for the inter ‐ planted area, with significant differences from the reference agricultural land (33.35 Mg ha − 1 ). Linear regression models were developed to predict the biomass and SOC at the tree scale, based on tree height ( R 2 >0.72 for biomass, and R 2 >0.62 for SOC). The overall on ‐ site C stock in the holm ‐ cork oak plantation was 35.11 Mg ha − 1 , representing a net C stock rise of 0.47 Mg ha − 1 yr − 1 . The ALS data allows a reliable estimation of C stocks in holm and cork oak plantations and high ‐ resolution maps of on ‐ site C stocks are useful for silvicultural planning. The cost of ALS data acquisition has decreased and this method can be generalised to plantations of other Mediterranean species established on agricultural lands at regional scales. However, an increase of filed data and the availability of local biomass and, in particular, SOC will improve accurate quantification of the C stocks from allometric equations, and extrapolation to large planted areas.


Introduction
Afforestation of set-aside agricultural lands was promoted by the European Economic Community 's approval of the CEE 2080/92 directive, which established a Community aid scheme for forestry measures in agriculture to promote the transformation of agricultural lands into forest [1]. The ecosystem importance of forests has long been recognised worldwide. From 1994 to 2006, an overall area of 232,858 ha of new forest plantations was established in Andalusia (Southwestern Spain), of which 83,000 ha were planted using holm oak (Quercus ilex L. subsp. ballota [Desf.] Samp.) as the main species [2]. Andalusia was the European region with the largest area of holm oak plantations in the Community forest plantation scheme due to the ecological and economic importance of this species in dehesa (savanna-type Mediterranean ecosystems). The Southwestern Iberian Peninsula is the area of maximum representation of such ecosystems, with 1,200,000 ha in Southern Spain.
Recent studies [3,4,5] indicated that agricultural abandonment may be an important and lowcost strategy for mitigation of climate change due to the vegetation recovery and increase in soil organic matter. Carbon is stored in forest ecosystems in different fractions, not only in the biomass of living vegetation but also in litter and forest soils [6]. Soils are the main terrestrial sink of atmospheric CO2, with storage in both organic and inorganic forms [6,7] and the largest proportion (above 50%) occurring in the upper horizons [8]. Soils store three-times more CO2 than living biomass and twotimes more than the atmosphere [9]. This C reservoir is not inert; it is immersed in the dynamics of the accumulation and decomposition of organic matter. About 70% of soil organic C (SOC) is stored under forests and the amount can be 1.5 times greater than that stored in the tree biomass [10]. The origins of C inputs into soils, particularly forest soils, are mainly related to plant litterfall and animal residua in different states of decomposition, parts of living or dead soil organisms, material synthesised by the soil microbial population and root residues (including exudates) [11]. The organic matter in the superficial layers of the soil has characteristics different from those of the organic matter in the deep layers. In the same way, the superficial and deep layers are not exposed to the same environmental factors. The SOC can be stored in soil for tens of centuries under suitable conditions, playing an essential role in the cycling of the nutrients required by plants and, therefore, in their growth potential [12]. However, the prevailing climatic conditions, land-use changes and land management, among other variables, can significantly affect the C storage capacity of a given ecosystem [13].
Land-use changes influences the C flows in natural and artificial ecosystems [5,12]. In this regard, afforestation of marginal agricultural lands is a relevant tool for the reversal of several soil degradation problems and improves C sequestration in plant biomass and soil [7], in spite of the fact that an increase in plant biomass production does not imply an increase in soil C stocks [13]. The assessment and quantification of the plantation C flux are essential to develop rapid, precise and economical measurements of soil C. Reliable quantification of the soil C storage and dynamic is also important to the identification of the most-appropriate policies and management practices for increasing, or at least stabilising, storage of SOC [14].
Several studies have described the relationship between SOC storage and the characteristics of the soils, climate, management and vegetation (tree species) [15,16]. Fernández-Getino et al. [17] emphasised the dependence of SOC on dasometric variables such as tree height, basal area and diameter in Quercus pyrenaica Willd. plantations. However, there is less information available on the relative importance of such dasometric variables to the C storage in different soil layers.
Recently, the assessment of forest resources and their dasometric characteristics has evolved significantly with the use of remote-sensing techniques, which have been replacing or are complementing traditional forest inventories. This is the case of airborne laser scanning (ALS), a remote-sensing technology based on LiDAR (laser imaging detection and ranging) that provides high-precision dasometric information on forest stands with considerable savings in inventory costs [18]. The ALS technique allows estimation of the main dasometric variables of a forest plantation (height, normal diameter, basimetric area, fraction of cover, density or volume of aerial biomass) from the generated point cloud, the calculation of its metrics and the application of algorithms that permit the development of models that relate these metrics to the main attributes of the forest [19].
In spite of the great advances in LiDAR technology and the precision that this technique offers in the estimations of the forest structure in forest inventories, research aimed at estimating the stock of SOC of Mediterranean plantations of Quercus spp. is scarce [20]. However, the estimation of this variable offers important information for the planning of the sustainable management of these forests to mitigate the effects of climate change and to increase the CO2 absorption capacity of Mediterranean forests. Additionally, the estimation of the contribution of afforestation to climate change mitigation requires knowledge on species-specific and country-specific biomass/C data and high-resolution regional models that contribute to improving the transparency and accuracy in the accounting for emissions/removals resulting from the land use, land-use change and forestry (LULUCF) sector, to fulfil the accounting obligations of European Union (EU) states [21].
In this study, high-density ALS data (12 points m -2 ) and different modelling techniques were used to estimate the height, aerial biomass and SOC of set-aside agricultural lands that have been afforested with a holm and cork oak (Quercus suber L.) plantation under Mediterranean conditions. The specific objectives were: i) to quantify the biomass and SOC present in a 14-year-old holm oakcork oak afforestation; ii) to develop allometric equations that relate the biomass and SOC stocks to tree height for the holm oak-cork oak plantation; iii) to estimate the on-site C stocks using highdensity ALS data and models based on regression techniques; and iv) to determine the total on-site C stocks in the holm oak-cork oak afforestation studied. The results of this study are discussed in the context of the use of ALS systems for the assessment of forest biomass and soil C stocks in the afforestation programmes of set-aside agricultural lands, with a view to evaluating the environmental services of climate change mitigation and CO2 fixation provided with an adequate accuracy and considerable savings in time and inventory costs.

Study Area
The study area is located in Puebla de Guzman (37°42'27" N-7°17'54" W, Huelva, Southwestern Spain, 190 m.a.s.l.). The terrain is undulating, annual rainfall averages 569 mm (2000-2010) and the annual mean temperature is 16.8 °C, according to the closest meteorological station (37°33ʹ07ʹʹ N-07°14ʹ54ʹʹ W, 4.95 km from the study site; Figure S1, Supplementary Material). Geologically, the area is associated with the Southern Portuguese zone, with abundant Devonian slates and shales. The soils of the study site are predominantly reddish-brown Eutric Regosols (Entisols) [22] with a low organic matter content (2.62%), low fertility and an almost-total absence of calcium carbonate. Soil pH (H2O 1:2.5) ranges from 5.6 to 5.8 and the texture from sandy loam to loam. The natural vegetation consists of a mixture of sparse evergreen shrubs (predominantly Cistus laurifolius L.) The holm-cork oak plantation was established in 1996 (it was 14 years old at the time of our measurements, 2010) in a 5 m × 6 m pattern, equivalent to a density of 330 trees ha -1 , although the density was reduced with time (to 208 trees ha -1 , 82% holm oak, 18% cork oak, Table 1), on 250 ha of degraded shrub pasture (Cistus ladanifer). This land was previously used for livestock production until it was afforested as part of the EU afforestation scheme. The planting was performed in two passes, using disc harrows drawn by a 70-hp farm tractor to remove shrub vegetation (20 cm depth) followed by linear subsoiling with a shank (40 cm depth) along the planting line. To control the spontaneous vegetation, tillage was performed every year. The plantations were never harvested, fertilised or irrigated. The trees were pruned when they were between 5 and 8 years of age, and the residues were chipped and dispersed in the ground.

Biomass Samples and Calculation
Our study used several data sets and required the development of remote-sensing indices and data analysis procedures. Therefore, a flow chart outlining the steps and relationships of each process is provided in Figure 1.
In October 2010, an extensive field survey was performed to identify a plot appropriate for the study ( Figure S1, Supplementary Material). The selected plot covers an area of 4900 m 2 (70 m × 70 m) of regular plantation, with an average slope of 20% and a NW exposure, in which all the selected trees were measured and soil samples were collected. The plot was selected based on the representative soil-forming factors for the study area corresponding to an average parental material, vegetation, and topographic conditions, and with sufficient internal variability (microtopoghaphy) to collect the condition variation of the whole plantation factors. The diameter at breast height (1.3 m above ground level -dbh-cm), diameter at the base of the trunk (Db, cm), crown diameter (Dc, m, as the average of the dimension of the canopy parallel to the plantation line and its perpendicular dimension and total height (distance between the root collar and the base of the terminal bud of the tree, H, m) of 308 trees (199 holm oak and 109 cork oak) were measured within the plot with a caliper (Haglöf, Sweden) and a Vertex III hypsometer (Haglöf, Sweden) ( Table 1). The trees were always selected in the central rows of the plantation, to avoid the edge effect.
During the forest inventory, data from the destructive sampling of 14 trees per species were used to develop biomass equations, based on the assumption that the allometry of stem biomass for trees growing under similar conditions does not differ significantly within the same genus [23]. These trees were uprooted and the biomass fractions were assessed in situ from direct measurements, by weighing on an industrial balance and, later, oven drying of biomass samples. At the tree base and at heights corresponding to 50% and 75% of the total height, two 2-cm-thick discs were extracted, which were weighed fresh (with bark) in the field and then dried to a constant weight at 103°C, to estimate the dry matter content (W, kg). Branch biomass was assessed on the basis of the fresh to dry weight ratio, which was estimated for branches selected from three segments of the tree crown that were later oven-dried. The sample branches from each tree were also used to estimate the foliage biomass. The aboveground biomass (Wa) was represented by the sum of the stem, branches and foliage. The total biomass (Wt) was represented by the sum of Wa and the belowground dry weight biomass (Wb). The biomass components were converted to C stocks using the 47.5% C fraction determined for Quercus by elemental analysis. A standard coefficient of 0.5 [23] was used to obtain the biomass C content.

Soil Carbon Stock Samples and Calculation
In May 2010, 45 soil profiles were taken systematically in the plot and in the reference agricultural land using an 8-cm-diameter steel auger. Within the plantation, 14 soil profiles were taken perpendicularly one metre away from the trunk (under crown sample), 28 were taken three metres away from the trunk (inter-planted sample) and three were taken in the reference agricultural land. Samples of each profile were extracted every 10 cm to a depth of 40 cm, because this is the general depth at which the parent material is found (N = 180 soil samples). The samples were bagged independently for transport to the laboratory, where they were processed according to the United States Department of Agriculture-Natural Resources Conservation System [24] methods. Briefly, the soil samples were air-dried and then rolled and sieved to separate the fine fraction (mesh size 2 mm). The resultant fine-earth fraction was stored for SOC determination. The gravel (>2 mm) was weighed and stored separately.
The organic carbon content of the fine-earth fraction was determined through wet oxidation by the Walkley and Black method [25]. The bulk density (BD) of the fine-earth fraction in every soil layer was estimated as follows (g cm -3 ) (Equation (1)): . .

,
( where W(fe) is the weight of the fine-earth fraction, Vol.cil. is the volume of the cylinder (10 cm in height), W(gr) is the weight of the gravel fraction and BD(gr) is the bulk density of the gravels, obtained in every case as the quotient between the weight of the gravels and its volume (obtained by the displacement of water).
The overall soil organic carbon stock (SOC40-S) in the first 40 cm from the soil surface was calculated by adding together the values obtained for each layer. A sub-metre global satellite receiver (Leica Zeno 20 GIS, Leica Geosystems, Switzerland) was used to survey plot centres and soil samples.
All data are represented as means ± standard error: differences between biomass stocks were tested by the Student´s t-test, and differences in SOC by one-way analysis of variance (ANOVA). A significant difference was accepted at P < 0.05. We performed all analyses using R software [28].

Tree Biomass and Soil Organic Carbon (SOC) Equations
All the variables were examined to ascertain whether they fitted a normal distribution (Shapiro-Wilk test, P < 0.05) and whether the variances were homogeneous (Levene test for the principal independent factor, P <0.05). When the data distribution did not fit a normal curve, the variable was subjected to a logarithmic transformation.
Biomass and SOC values were regressed against individual dasometric variables to obtain Pearson´s coefficients. However, to estimate total biomass per tree (Wi) and SOC-S, only total H was used as a predictor variable to obtain allometric coefficients; thus, it was possible to make a fast estimation of the C stocks by means of ALS data. Linear, exponential, power and logarithmic regressions models were used. The Akaike information criterion (AIC) was used for the comparison of models [29], choosing the equation with the lowest ΔAIC statistic (the difference between AIC values for two nested models), as this can be interpreted as a trade-off between model fit and complexity. The coefficient of determination (R 2 ) and root mean square error (RMSE) were considered as secondary criteria of model evaluation. The R 2 value was used to assess model fit for equations with one parameter, while equations with multiple parameters were assessed using the adjusted coefficient of determination (R 2 adj).
All the models reported were tested using the Global Validation of Linear Models Assumptions. package of R software [30], which performs a global validation of linear model assumptions.

Aerial Laser Scanning (ALS) Data and Processing
The ALS data were acquired on 6 April 2013 by the company Heliografics Fotogrametria S.L. (Alicante, Spain), using an ALS60 laser scanner (Leica-Geosystems AG, Heerbrugg, Switzerland). The scan frequency was 100 Hz, the laser pulse repetition rate was 187.6 kHz, the illuminated footprint diameter was 28 cm and the field-of-view (FOV) was 13° (Table S2, Supplementary Material). The field was scanned by a plane from a flight altitude of 1300 m. The ALS data were acquired with a point density of 12 points m -2 . They were geo-referenced in the European Terrestrial Reference System 1989 (ETRS89) coordinate system. The planimetric coordinates (x and y) and ellipsoidal height values were computed for all echoes. The 3-year delay between the ALS data acquisition (2013) and the field data collection (2010) was not considered a significant source of error, as the plantations under study did not change considerably during that period (mean diametric annual increment = 3.46 ± 0.03 mm yr -1 for holm oak and 4.92 ± 0.04 mm yr -1 for cork oak) [31]. The forest-stand homogeneity and geographic distribution make this statistic robust and informative.
A pit-free canopy height model (CHM) was generated by following the steps in Khosravipour et al. [32], using Lastools [33]. Then, a region growing segmentation algorithm [34] was applied for the delineation of tree crowns. This algorithm requires a seed point input to locate the treetops, which are retrieved with local maxima-filtering approaches. Once the seeds are located, neighbouring pixels are compared with the initial point and merged until some specified threshold criterion is reached. The input parameters of variance in feature space, variance in position space and similarity threshold were tested, evaluated and selected by following the methodology in Varo-Martínez et al. [34]. For this study, processing was undertaken in the SAGA (System for Automated Geoscientific Analyses) geographic information system (GIS) [35].
The ALS-based height metrics obtained in each plot were: minimum, maximum, mean, median, standard deviation, variance, coefficient of variation, interquartile distance, skewness, kurtosis, AAD (average absolute deviation), L-Moments (1-4) and percentile values (H_P5 to H_P95 in five-unit intervals and H_P99) [36] (Table S2, Supplementary Material). For further details of the procedure used to obtain such ALS metrics, see the steps described in González-Ferreiro et al. [37]. The metrics involved in the equation of the H model were obtained using the "GridMetrics" and "CSV2Grid" commands implemented in FUSION LDV 3.30 [38]. All ALS indices were calculated with lidR package v2.0.0 [39].

Laser Imaging Detection and Ranging (LiDAR) Analysis and Height Models
Relationships among ALS metrics and H were determined using stepwise regression models. Explanatory variables were selected, firstly using variance inflation factor (VIF) analysis to check multicollinearity and secondly choosing the variable or combination of variables which minimised the generalized root mean square distance when variables were added or deleted one at a time [40]. Thus, once the best predictor variables had been selected, we used stepwise regression to find the best variables for H estimation [37]. In total, 205 segments (trees) (2/3 of the data) were assigned to the building of a predictive classification model using ALS metrics and tree height. We used the same criteria of model evaluation used for the tree and SOC models. We performed all analyses with R software, version 3.4.0, using the yaImpute R package [40] to calculate variable selection and regression and the usdm package to perform collinearity analysis [41].

Cartography of C Stocks
Finally, a C stock map of the holm-cork oak cover in the study area was generated. The models obtained in Section 2.5 were used to calculate H from the ALS metrics at the tree scale (250 ha, 51,093 trees). Once the H of each tree had been calculated, the allometric models considering H as an independent variable were applied to determine Wt-S and SOC40-S.
Maps were generated at the tree scale to allow estimation of the on-site C stocks per hectare for each species and were summarised in a tabular format for height classes. Ratios of biomass and SOC C stocks were calculated for the 14-year time period. With regard to CO2 absorption and economic balance, we considered the current price on the C market (24.77 € (Mg CO2) -1 .

C Stock in Biomass after Afforestation
The biomass C stocks for all fractions on the afforested land are presented in Table 1. The contributions of different components to the total tree biomass varied considerably between holm and cork oak and their fractions, over the study period (14 yr). The C stock in biomass (Wt-S) was lower in holm oak (27.9 kg C tree -1 , P < 0.001) than in cork oak (41.1 kg C tree -1 ). The Wa-S accounted for most of the total tree biomass, contributing 54.1% and 75.4%, respectively, for holm and cork oak, while the respective root:shoot ratios were 0.84 and 0.35.

SOC Stock after Afforestation
The SOCt-S under the agricultural land was 33.35 Mg ha -1 , similar to that under the tree crown (35.60 Mg ha -1 ) and significantly (F = 5.08, P = 0.011) higher than in the inter-plantation area (25.34 Mg ha -1 ) ( Table 1). Additionally, afforestation gave rise to variation in the proportions of SOC-S present along the soil profile, particularly in the deeper layers (Table 1, Figure S2, Supplementary Material).

Individual Biomass and SOC Equations
Table S2 (Supplementary Material) shows the Pearson's correlation coefficients (r) that were obtained between biomass and soil C stocks and dasometric variables. Total height was the best individual predictor of biomass (Wt-S, r = 0.85), for holm-cork oak, and SOC (SOC40-S, r = 0.85). Other dasometric variables also had correlation coefficients that showed relationships with biomass and SOC that were significant, but less so.

Prediction of Height from LiDAR Data
Following the selection of the independent variable data, the H model used the height (H_P90 = Elevation percentile 90; H_P60 = Elevation percentile 60; Elev.sqrt = Elevation quadratic mean; Elev.mean = Elevation mean; Elev.MAD = the median of median absolute deviation; Elev.AAD = average absolute deviation), the returns (bin70 = Fraction of return points between the 70 th percentile height and the maximum height) and the canopy variables (Can.rel = Canopy relief ratio) ( Table 4,  Table S2, Supplementary Material). The predictor variables had a high correlation with H, which suggests that LiDAR metrics are a good representation of the forest metrics. Table 4 shows the best H-estimation models, based on the stepwise regression method. Overall, the agreement for the holmcork oak trees was good, the bias not differing significantly from zero (t-test, P < 0.05). The RMSE was small (0.58 m, 3.98%) considering that these approaches were compared at the tree level, where the inter-tree variability is high.

Total On-Site C Stock in the Plantation
Using the segmentation process, a total of 51,093 trees were delineated, with an average height of 2.53 m ( Table 5). The best regression models were selected to spatially estimate the C stocks for the study area (Wt-S + SOC40-S, see Table 2 and Table 3, Figure 2 and Figure 3). Based on the tree cartography and these models, the distribution maps with predicted values of C stocks in the planted area show a mosaic of C stock patterns in the holm-cork oak plantation (Figure 2 and Figure 3). We obtained a mean Wt-S value of 4.89 Mg ha -1 , and a mean SOC40-S value of 25.31 Mg ha -1-ranging from 34.17 Mg ha -1 , under the tree crown, to 24.94 Mg ha -1 , in the inter-plantation area (Table 5). In the holm-cork oak plantation, the on-site C stock (Wt-S + SOC40-S) reached 35.11 Mg ha -1 , which translates into a landscape-wide mean annual rate of C sequestration of 0.47 Mg ha -1 yr -1 for the plantation (Table 5).

Effects of Afforestation on Biomass Carbon and Soil Organic Carbon Stocks
The average overall biomass C stock calculated was 4.89 Mg ha -1 in the holm-cork oak plantation. Our values are higher than those reported for low-density plantations with the same species (1.14 ± 0.13 Mg ha -1 , Q. ilex; 0.85 ± 0.11 Mg ha -1 , Q. suber; [42]) but lower than those reported for similar species (6.78 ± 2.68 Mg ha -1 , Ceratonia siliqua; [42]) at the same age (12-14 yr). These discrepancies can be explained by the differences in plantation density (100-300 trees ha -1 ) [43]. Aboveground biomass accounted for the largest proportion of individual tree biomass, for both species. This division agrees with results from studies on the growth of Quercus species in Spain [44], namely an aboveground biomass proportion greater than 60%. The lower investment in belowground biomass of the trees in the surveyed plantation indicates soil limitations or a root morphological pattern related to seedling quality, sacrificing root growth for growth in height and an expanding crown area [45]. However, the belowground biomass is still an important C pool for the studied species, with a relatively-high proportion of the total tree biomass determined in this study.
Regarding the SOC-S, in our study SOC40-S was 33.35 Mg ha -1 and 36.90 Mg ha -1 on agricultural land and under the tree crown, respectively, significantly higher than in the inter-plantation area (29.26 Mg ha -1 ). The SOC stock decreased along the soil profile, which coincides with the normal distribution in this type of plantation [6]. Our values are higher than those reported in other studies, including a study of mature semi-natural populations of Quercus suber, in which the C content of the first 50 cm of soil was calculated to be 24.20 Mg ha -1 [46], and a study of plantations of a similar species (Quercus pyrenaica), which reported a value of 33 Mg ha -1 [17,47]. The differences among species may be due to differences in litter production, litter quality and microclimate-which may result in different proportions of leaf organic C in the increased SOC [48].
In the literature, most studies reported an increase in SOC-S following reforestation or afforestation of cultivated land [49,50]. In contrast, other studies reported decreases or non-significant changes in SOC-S after conversion from agriculture to forest plantation (e.g., Vesterdal et al. [51]). As a consequence of the 14-year-old plantation, the SOC was 10.64% higher under the tree crowns than in the agricultural land, but 12.26% lower in the inter-plantation area than in the agricultural land. The positive difference in SOC-S between the soil under the tree crown and that of the agricultural land, in particular in the shallowest soil layer (0-20 cm), may be related to several factors such as i) variation in the quality and quantity of inputs (litterfall and logging residues), ii) greater density of roots, iii) an increase in N use efficiency, iv) greater diversity and abundance of soil organisms (bacteria, fungi, earthworms and others), v) protection of both new and old SOC in organic-mineral complexes within aggregates (<50 μm) and vi) a greater micro-climatic effect on litter decomposition, resulting in rapid SOC accumulation and long-term storage in the soil [52][53][54]. These findings support the idea that holm-cork oak plantations on former agricultural land can accumulate large amounts of SOC and reach an equilibrium during the first two decades after their establishment [6]. However, the observed losses of SOC in the continuously-cultivated inter-plantation area, compared to the agricultural field, could be due to several factors: (i) different C sources in the agricultural crops (i.e., a higher root biomass contribution, [55]), (ii) increased decomposition rates in cultivated soils due to higher soil temperatures, (iii) the different soil management practices (e.g., mechanical weed control in the holm-cork oak plantation, [56]), (iv) the time elapsed since the establishment of the plantation and (v) increased infiltration and, thereby, increased loss of SOC such as DOC (dissolved organic carbon), because of reduced rainfall interception [50].

Individual Biomass Equations
The allometric equations obtained in this study permit the biomass and soil C stocks to be predicted as a function of tree height. Total height alone was the independent variable that best described the different biomass components evaluated for a holm-cork oak plantation (R 2 > 0.79; p-DW and p-WT <0.05). Belowground biomass showed a similar fit, not confirming previous reports showing that Wb is a major component of uncertainty when measuring total tree biomass [57]. Similarly-good fits of high R 2 have been reported by Ruiz Peinado et al. [44] for several forest tree species in Spain, using height and other predictor variables (e.g., dbh). The relationship was less pronounced for overall biomass (R 2 > 0.62; p-DW and p-WT < 0.05), which makes sense considering that tree form and pruning modify the canopy size of individual trees and alter more the allometric relationships [58]. Also, this variability may be related to the difficulty in accurately estimating the biomass of small trees. Despite these limitations, H was a useful parameter to estimate the biomass C stock in low-density Quercus plantations. Other variables, such as dbh or crown diameter, were not included as additional predictors of tree biomass components in this study because we were looking to relate tree biomass to ALS data [59].

SOC Equations
Regarding soil C stocks, an equation was generated to predict the biomass and SOC40 C stocks based on tree height (R 2 = 0.69 for holm and cork oak together). Several authors have related the soil C stock to the age of the plantation and therefore to its growth [60]. However, to the best of our knowledge, no study has used the height to estimate the SOC for holm and cork oak. This dependence may arise from the fact that the main source of incorporation of organic matter into the soil is through litter fall [61]. Additionally, due to the pruning treatments practiced in Quercus plantations, soil C stocks are strongly associated with height-as in the case of other species that undergo regular management [9], associated with pruning and thinning [62].

Prediction of Height from ALS Data
In recent years, ALS sensors have been used mainly to estimate forest parameters, being a critical technology for forest inventories and management in natural [63] and planted forests [20]. According to the regression models generated in this study, height metrics (H_P60, H_P90) and another parameter concerning the horizontal distribution of the point cloud (i.e., the canopy density or the percentage of returns above a certain height threshold) can be used to generate accurate information on tree height in fairly homogeneous holm-cork oak plantations. Several authors have estimated the height in Pinus plantations using ALS data [37,64], and the height in Q. suber plantations in Portugal has been estimated using unmanned aerial vehicle data based on similar ALS metrics [65]. The coefficients of determination between the ALS-derived metrics and tree height obtained here are comparable to those of other studies in forest plantations (from 0.60 to 0.90; [36,66,67]). The main sources of error in our study may be related to the time delay of the ALS data, the precision in the geolocation of individual trees and the tree allometry. Working with LiDAR data demands high precision in field plots coordinates, for the comparison of dasometric data with ALS metrics derived from the point cloud. Thus, higher values of RMSE can be due to problems with individual tree geolocation (sometimes with errors greater than 1 m), such that the LiDAR point clouds do not correspond exactly with the dasometric values obtained through the field inventories. Finally, a high degree of tree overlapping and a tendency to produce many shoots in young (small) Quercus trees may result in detection errors at all heights [65].
The high cost of fieldwork is usually prohibitive in the production of detailed tree-height maps for large plantations. The ALS sampling methods used across planted areas have improved our ability to monitor forest parameters, biomass and C content. Low-density ALS data [20,37], sometimes with less than 1 point m -2 , allow the use of ALS techniques in the development of regionalscale C inventories based on tree height information at the tree and stand scales [20,68]. However, because the density of the ALS data is potentially the main factor limiting tree height estimation in forest plantations in agricultural lands [69], and we have used high-density ALS data (12 points m -2 ), further investigation is required to develop systematic monitoring of canopy density and to relate this variable to ALS height estimations.

Total On-Site C Stock and Cartography
Finally, we created a set of C stock maps of the holm-cork oak plantation at the tree scale (51,093 trees) to estimate the on-site C stock. By applying the allometric models generated in this study, the biomass and soil C stocks were evaluated indirectly and efficiently for the study area. Our estimate of Wt-S + SOC40-S (35.11 Mg ha -1 , Table 4) shows low growth rates and low SOC, and thus small C stocks, in comparison with other plantations on former agricultural land [63]. However, when the net balance of C absorption is calculated, it is evident that the holm-cork oak plantation was more efficient than the reference agricultural land, especially under the tree crown, where the Wt-S + SOC40-S sequestration was 41.79 Mg ha -1 ( Table 1). The overall on-site C stock (Wt + SOC40-S) was 8777.50 Mg in the 14-year-old, 250-ha plantation. When considering this overall C stock, we obtained a positive C accumulation rate after establishment of 0.475 Mg ha -1 yr -1 with respect to the agricultural land. Segura et al. [70] found on-site C stock changes in the SOC40 of 0.61 Mg ha -1 yr -1 during the first two decades following the planting of Pinus halepensis, although this was in a semi-arid Mediterranean area. The positive shift in the on-site C stock in our study indicates that changes in C pools after tree planting were mainly due to increases in tree biomass and in litter deposition under the tree crown, resulting in a net accumulation rate of 0.34 Mg ha -1 yr -1 of biomass and 0.12 Mg ha -1 yr -1 of SOC40 (Table 4). However, these inputs of C stocks did not compensate the losses of SOC in the inter-plantation area. Afforestation of agricultural lands may result in faster on-site C stocks accumulation if adequate soil and vegetation management are applied, suggesting that low-density forest plantations have high potential for C sequestration in Mediterranean areas. To support this forest management plan, our results show that accurate mapping of the C stock distribution can be modelled with good precision using ALS data and a limited number of field measurements, as observed in previous studies [20].
Our results outline the important role that afforestation on agricultural land may play in increasing the C sink in the coming years. Under a non-intervention scenario (agricultural land), SOC C sequestration would be 30,348.50 Mg CO2 (for an equivalent area, 250 ha), only slightly lower than in the afforestation scenario (14-year old plantation, 31,958.18 Mg CO2, an increment of 5.30%) due to the SOC losses in the inter-plantation area. Thus, given the potential contribution of afforestation to C sequestration under adequate management practices (e.g., soil conservation practices and vegetation management), its implementation must be encouraged to counter atmospheric CO2 emissions [71]. New afforestation efforts on marginal agricultural lands should be made, following schemes similar to those of the EU. This target should be included in the rural development strategy of the EU, thus ensuring the maintenance of ecosystem services related to planted forests. Many previous plantations, established in the EU afforestation scheme, have been abandoned and future management should prioritise on-site C sequestration in long-lived pools, as well as conservation and extensive soil-vegetation management practices to increase potential C sequestration [72].
One way to estimate the economic value underlying C sequestration is based on the CO2 stock exchange in this potential scenario. In the most-optimistic situation (an accumulation rate of 0.475 Mg ha -1 yr -1 ), and considering the current holm-cork oak planted area in Andalusia (37,670 ha) with an average age of 14 years (1996-2010), a total of 65,131.43 Mg CO2 yr -1 could be sequestered, almost 0.12% of the current fossil fuel emissions of Andalusia (equivalent to 52,310 kMg CO2). The economic value of C sequestration would, therefore, be 1,613,305.52 € yr -1 in the intervention scenario (24.77 € (Mg CO2) -1 in the EU Emission Trading averaged monthly value for the year 2019). In Andalusia, 113,952 ha of forest have been planted in the EU scheme, and represent an alternative to increase the capacity of this area to provide ecosystem and economic services, goods and woodland conservation if these benefits are returned to the owners as subsidies and incentives for forest plantation management [73,74].

Conclusions
ALS is one of the most powerful tools to estimate C-stock in forested long areas. ALS offers precise data collection methods allowing data to be collected faster, with greater precision and at lower cost than field methods. On the contrary, ALS data are expensive, have a high point density and are difficult to process for large volumes of data, but low-density national ALS covers can help to minimise this limitation. Finally, the use of relationships between the ALS data and C biomass and soil require equations that must be parameterized locally. However, on this, and other previous studies, ALS has proven to be a very effective way to estimate C-stock in Mediterranean Quercus plantations under the European Community aid scheme for forestry measures in agriculture. The change from agricultural land to a forest plantation resulted in an increase in the biomass C stock. However, it also resulted in a slight loss of SOC 14 years after afforestation. Afforestation of agricultural lands may result in positive on-site C stocks accumulation rates if adequate soil and vegetation management are applied, especially in the inter-plantation area. Adjusted allometric models have provided accurate estimation of biomass and SOC under local conditions. Based on these models, this paper proposes a methodological approach to use ALS data as an indirect estimator of tree biomass and SOC in forest plantations in agricultural lands. The approach relies on the estimation of canopy height using ALS metrics at the tree scale. The detailed cartography obtained with ALS data represents a reliable picture of the current situation of C stock pools in holm and cork oak plantations. These high-resolution maps of on-site C stocks are useful for silvicultural planning and future monitoring of potential changes in this area, in terms of C sequestration in forests, due to sustainable silvicultural practices. Recently, the cost of ALS data acquisition has decreased-they can be obtained for free in some countries (e.g., in Spain, 0.5-1.0 points m −2 )-and is comparable to or even less than the cost of large-scale field data collection. This method can be generalised to plantations of other Mediterranean species established on agricultural lands at regional scales, reducing the costs of C stocks inventories. Our results also provide reference values for the species planted most in EU-afforestation agricultural land schemes in Southern Spain.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1: Figure S1. Location of the holm-cork oak (Quercus ilex and Q. suber) plantation in Puebla de Guzmán (Southwestern Spain) and climatic characterisation of the study area. Figure S2. Soil organic carbon content along soil profiles in the holm-cork oak plantation. Table S1. Pearson's coefficients (r) between the weighted SOC (kg m -2 ); and the dasometric values of total height (ht), canopy height (hc), normal diameter at breast height (dbh), base diameter (db) and canopy diameter (dc). Table S2. Correlation coefficients describing the strength of the linear relationships between the stand metrics (mean height, mean diameter, basal area, stand density and biomass) obtained with laser imaging detection and ranging (LiDAR) and the measured stand metrics.

Acknowledgments:
The authors thank the Andalucía Department of Agriculture and Environment, which provided access to and background information on the field site, in particular Fernando Piñón. We are very grateful to Rafael Sánchez de la Cuesta, Francisco Pereña, Emidio Silveiro, and Andrés Cortés for their valuable assistance during field work and data acquisition and processing; without them, this project would not have been possible. We also thank Marta Álvarez for guidance regarding interpretation of the soil data. We also acknowledge the institutional support of the University of Cordoba-Campus de Excelencia CeiA3.

Conflicts of Interest:
The authors declare no conflict of interest.